Commit 85ba62d4 authored by roelo008's avatar roelo008

Berekent aantal inwoners onder/boven de norm in situatie zonder vegetatie

parent 4c3e413b
'''
"""
ESD Model 2J - fijnstof-concentratie in de lucht.
code: Hans Roelofsen
idee: Marjolein Lof
......@@ -29,7 +29,7 @@ conc_scen = (1-0.137 * myear) * conc_blanco (note output <= input)
Aannames
input rasters zijn geharmoniseerd, ie identieke Affine
'''
"""
import os
import datetime
......@@ -38,12 +38,18 @@ import rasterio as rio
from z_utils import esd
'''INPUT'''
# input variables can be changed to command line arguments later on!
'''INPUTS. Can be changed to command line arguments later on!'''
# Rasters
conc_rivm_in = r'c:\apps\proj_code\2J\input\rpm25_2016_resamp.tif' # path to RIVM referentie concentratie
lceu_orig_in = r'x:\ESD_SpatialData\EcosysteemUnits_LCEU\f_non-BIG-tiff\a_orig\LCEU10_orig_gdal.tif' # orig LCEU map
lceu_scen_in = r'x:\ESD_SpatialData\EcosysteemUnits_LCEU\f_non-BIG-tiff\e_hdbb\LCEU10_hdbb_gdal.tif' # LCEU scenario
cbs_inwoners = r'x:\ESD_SpatialData\CBS_4kant_gegevens\rasterized\CBS_VK100_2018_v1.tif' # path to CBS inwoners grid
# Output dir
out_base_dir = r'c:\apps\proj_code\esd\c_out'
# Parameters
scen_naam = 'hdbb' # scenario naam
norm = 10 # norm voor concentratie
px = 10 # pixel grootte in m van de LCEU kaart en eigenlijk ook van alle andere rasters.
......@@ -59,12 +65,12 @@ cbs_bevol = rio.open(cbs_inwoners)
conc_rivm_arr = conc_rivm.read(1) # numpy array
lceu_orig_arr = lceu_orig.read(1)
lceu_scen_arr = lceu_scen.read(1)
# div by 100 to compensate resampling from 100m to 10m pixels. Should actually be done in QGIS preprocessing
cbs_bevol_arr = np.divide(cbs_bevol.read(1), 100).astype(np.float32)
'''Check Affine consistency'''
if esd.check_affine_consistency(conc_rivm, lceu_orig, lceu_scen, cbs_bevol):
pass
esd.check_affine_consistency(conc_rivm, lceu_orig, lceu_scen, cbs_bevol)
'''Calculate blanco concentratie '''
myear_blanco = esd.vec_translate(a=lceu_orig_arr, my_dict=myear)
......@@ -74,31 +80,38 @@ conc_blanco = esd.conc_blanco(conc_rivm_arr, myear_blanco)
myear_scen = esd.vec_translate(a=lceu_scen_arr, my_dict=myear)
conc_scen = esd.conc_scenario(conc_blanco, myear_scen)
'''Treshold conc_scen to the norm'''
'''Treshold blanco concentratie andd scenario concentratie to the norm'''
eval_blanco = np.where(conc_blanco <= norm, 1, 0).astype(np.uint8)
eval_scen = np.where(conc_scen <= norm, 1, 0).astype(np.uint8) # set to 1 where value <= norm
'''Calculate population living below/above norm'''
# Done on the fly later on
'''Output numbers '''
true_pxls = eval_scen.sum() # sum, i.e. pixel count where value <= norm
true_area = np.product([true_pxls, np.divide(np.square(px), 1000000)]) # multiply pxl count with area per pixel
true_perc = np.product([np.divide(true_area, nl_onshore), 100]).round(2) # divide by reference area
# totale populatie in NL
total_pop = cbs_bevol_arr.sum().astype(np.uint64)
below_pop = np.where(eval_scen == 1, cbs_bevol_arr, 0).sum().astype(np.uint64)
above_pop = np.subtract(total_pop, below_pop)
below_perc = np.product([np.divide(below_pop, total_pop), 100]).round(2)
above_perc = np.product([np.divide(above_pop, total_pop), 100]).round(2)
# scenario population below and above treshold as number and percentage
scenario_below_pop = np.where(eval_scen == 1, cbs_bevol_arr, 0).sum().astype(np.uint64)
scenario_above_pop = np.subtract(total_pop, scenario_below_pop)
scenario_below_perc = np.product([np.divide(scenario_below_pop, total_pop), 100]).round(2)
scenario_above_perc = np.product([np.divide(scenario_above_pop, total_pop), 100]).round(2)
# w/o vegetation population below and above treshold as number and percentage
blanco_below_pop = np.where(eval_blanco == 1, cbs_bevol_arr, 0).sum().astype(np.uint64)
blanco_above_pop = np.subtract(total_pop, blanco_below_pop)
blanco_below_perc = np.product([np.divide(blanco_below_pop, total_pop), 100]).round(2)
blanco_above_perc = np.product([np.divide(blanco_above_pop, total_pop), 100]).round(2)
'''Write outputs to file'''
timestamp_brief = datetime.datetime.now().strftime("%Y%m%d%H%M")
timestamp_full = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
out_base_dir = r'c:\apps\proj_code\esd\c_out'
out_dir = os.path.join(out_base_dir, '{0}{1}'.format(timestamp_brief, scen_naam))
os.mkdir(out_dir)
out_dir = os.path.join(out_base_dir, '2J', '{0}{1}'.format(timestamp_brief, scen_naam))
os.makedirs(out_dir)
basename = 'ESD2J_rpm25_{0}'.format(scen_naam)
# Float64 rasters
prof = lceu_orig.profile
prof.update(dtype=rio.float64)
with rio.open(os.path.join(out_dir, basename+'_conc.tif'), 'w', **prof) as dst:
......@@ -107,13 +120,21 @@ with rio.open(os.path.join(out_dir, basename+'_conc.tif'), 'w', **prof) as dst:
with rio.open(os.path.join(out_dir, basename+'_conc_wo_veg.tif'), 'w', **prof) as dst:
dst.write(conc_blanco, 1) # write concentration w/o vegetation
# Float 32 rasters
prof.update(dtype=rio.float32)
with rio.open(os.path.join(out_dir, basename+'_pop_below_tresh.tif'), 'w', **prof) as dst:
with rio.open(os.path.join(out_dir, '{0}_scenario_pop_below_tresh.tif'.format(basename)), 'w', **prof) as dst:
dst.write(np.where(eval_scen == 1, cbs_bevol_arr, 0), 1) # write populating living below treshold raster
with rio.open(os.path.join(out_dir, basename+'_pop_above_tresh.tif'), 'w', **prof) as dst:
with rio.open(os.path.join(out_dir, '{0}_scenario_pop_above_tresh.tif'.format(basename)), 'w', **prof) as dst:
dst.write(np.where(eval_scen == 0, cbs_bevol_arr, 0), 1) # write population living above treshold raster
with rio.open(os.path.join(out_dir, '{0}_blanco_pop_below_treshold.tif'.format(basename)), 'w', **prof) as dst:
dst.write(np.where(eval_blanco == 1, cbs_bevol_arr, 0), 1)
with rio.open(os.path.join(out_dir, '{0}_blanco_pop_above_treshold.tif'.format(basename)), 'w', **prof) as dst:
dst.write(np.where(eval_blanco == 0, cbs_bevol_arr, 0), 1)
# Uint8 rasters
prof.update(dtype=rio.uint8)
with rio.open(os.path.join(out_dir, basename+'_below_tresh.tif'), 'w', **prof) as dst:
dst.write(eval_scen, 1) # write binary above/below treshold rasters
......@@ -125,16 +146,22 @@ footer = '\nMade with Python 3.5 using pandas, geopandas, by Hans Roelofsen, WEn
msg1 = 'For scenario {0} (based on {1}) there are {2} pixels w value <= {3}, which equals to:' \
'{4} square km, which is {5}% of onshore area of the Netherlands ({6} km2).\n'\
.format(scen_naam, lceu_scen_in, true_pxls, norm, true_area, true_perc, nl_onshore)
msg2 = 'There are {0} people living in the Netherlands in 2018 according to CBS, of which {1} ({2}%) live in areas ' \
'where the treshold is exceeded and {3} ({4}%) in areas below the treshold.' \
.format(total_pop, above_pop, above_perc, below_pop, below_perc)
msg2 = 'Anno 2020 there are {0} people living in the Netherlands in 2018 according to CBS. In the scenario ' \
'{1} ({2}%) would live in areas where the treshold is exceeded and {3} ({4}%) in areas below the treshold.\n' \
.format(total_pop, scenario_above_pop, scenario_above_perc, scenario_below_pop, scenario_below_perc)
msg3 = 'In a situation without vegetation, {0} people ({1}%) would live in areas where the treshold is exceeded ' \
'and {2} ({3}%) in areas below the treshold.\n'.format(blanco_above_pop, blanco_above_perc, blanco_below_pop,
blanco_below_perc)
print(msg1)
print(msg2)
print(msg3)
with open(os.path.join(out_dir, '{0}_report_{1}.txt'.format(os.path.splitext(scen_naam)[0], timestamp_brief)), 'w') \
as f:
f.write(msg1)
f.write(msg2)
f.write(msg3)
f.write(footer)
......
......@@ -17,7 +17,7 @@ def myear_tab():
def conc_blanco(ref, myear):
'''
"""
Bereken concentratie fijnstof zoals die zou zijn ZONDER vegetatie, gegeven referentie-concentratie
ref: array met referentie-concentratie, e.g. RIVM 2016
myear: array met afvangcapaciteit per pixel, e.g. combinatie LCEU, myear_Tab
......@@ -26,7 +26,7 @@ def conc_blanco(ref, myear):
out = ref * (1 + noemer/deler)
out = ref * (1 + foo)
out = ref * bar
'''
"""
noemer = np.multiply(myear, 0.137) # array
deler = np.subtract(1, noemer) # array
......@@ -37,28 +37,32 @@ def conc_blanco(ref, myear):
def conc_scenario(blanco, myear):
'''
"""
Bereken concentratie fijnstof gegeven een blanco situatie en de afvangcapaciteit zoals die is in een scenario
blanco: array met blanco concentratie
myear: array met afvangcapaciteit per pixel, e.g. combinatie LCEU, myear_Tab
conc_scen = (1-0.137 * myear) * blanco
'''
"""
out = np.multiply(np.subtract(1, np.multiply(0.137, myear)), blanco)
return out
def vec_translate(a, my_dict):
'''
"""
Function to apply my_dict to each element in in array a. See:
https://stackoverflow.com/questions/16992713/translate-every-element-in-numpy-array-according-to-key
'''
Note that key errors are not caught!
"""
return np.vectorize(my_dict.__getitem__)(a)
def check_affine_consistency(*rios):
'''Check if pixel size, origin, ncol and nrow for >=2 rasterio objects are identical. Return true or False + msg'''
"""
Check if pixel size, origin, ncol and nrow for >=2 rasterio objects are identical. Return true or False + msg
"""
# unpack arguments
rio_lst = [x for x in rios]
......@@ -74,7 +78,7 @@ def check_affine_consistency(*rios):
return True
else:
rasts = [os.path.basename(x.name) for x in rio_lst]
print('Inconsistent extent, pixel size and/or shape for provided rasters!\n\n')
print('\n'.join(['{0}: {1}'.format(a, b).replace('\n', '') for (a, b) in zip(rasts, affs)]))
return False
msg1 = 'Inconsistent extent, pixel size and/or shape for provided rasters!\n\n'
msg2 = '\n'.join(['{0}: {1}'.format(a, b).replace('\n', '') for (a, b) in zip(rasts, affs)])
return msg1, msg2
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment