...
 
Commits (2)
......@@ -2,6 +2,7 @@
ESD Model 2J - fijnstof-concentratie in de lucht.
code: Hans Roelofsen
idee: Marjolein Lof
datum: december 2019
Berekent concentratie fijnstof in de lucht gegeven landgebruik volgens een bepaald scenario
......@@ -38,9 +39,8 @@ from z_utils import esd
# input variables can be changed to command line arguments later on!
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' # path to original LCEU map
lceu_scen_in = r'x:\ESD_SpatialData\EcosysteemUnits_LCEU\f_non-BIG-tiff\b_huidig\LCEU10_huidig_gdal.tif' # path to LCEU scenario
# TODO: LCEU Huidig KeyError 128!
scen_naam = 'huidig' # scenario naam
lceu_scen_in = r'x:\ESD_SpatialData\EcosysteemUnits_LCEU\f_non-BIG-tiff\e_hdbb\LCEU10_hdbb_gdal.tif' # path to LCEU scenario
scen_naam = 'hdbb' # scenario naam
norm = 10 # norm voor concentratie
px = 10 # pixel grootte in m
myear = esd.myear_tab() # invangcapaciteit per landgebruiks-categorie
......@@ -64,7 +64,7 @@ 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'''
eval_scen = np.where(conc_scen <= norm, 1, 0) # set to 1 where value <= norm
eval_scen = np.where(conc_scen <= norm, 1, 0).astype(np.uint8) # set to 1 where value <= norm
'''Output numbers '''
true_pxls = eval_scen.sum() # sum, i.e. pixel count where value <= norm
......@@ -72,13 +72,24 @@ true_area = np.product([true_pxls, np.divide(np.square(px), 1000000)]) # multip
true_perc = np.product([np.divide(true_area, nl_onshore), 100]) # divide by reference area
'''Write outputs to file'''
# TODO
prof = lceu_orig.profile
out_dir = r'.\c_out'
basename = 'ESD2J_rpm25_{0}'.format(scen_naam)
prof.update(dtype=rio.float64)
with rio.open(os.path.join(out_dir, basename+'_conc.tif'), 'w', **prof) as dst:
dst.write(conc_scen, 1) # write concentration
prof.update(dtype=rio.uint8)
with rio.open(os.path.join(out_dir, basename+'_tresh.tif'), 'w', **prof) as dst:
dst.write(eval_scen, 1)
'''Write summary as message and to file'''
msg = '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)'\
.format(scen_naam, os.path.basename(lceu_scen_in), true_pxls, norm, true_area, true_perc, nl_onshore)
print(msg)
with open(r'.\{0}_report.txt'.format(os.path.splitext(scen_naam)[0]), 'w') as f:
with open(r'.\c_out\{0}_report.txt'.format(os.path.splitext(scen_naam)[0]), 'w') as f:
f.write(msg)
......
......@@ -3,11 +3,14 @@
import numpy as np
def myear_tab():
'''
For use in model 2J Luchtfiltering
:return: dictionary mapping relation b'teen LCEU land cover category and fine-dust particle cleasing capacity
'''
lceu = [0, 1, 2, 3, 4, 5, 6, 11, 12, 13, 21, 22, 23, 24, 25, 26, 27, 28, 29, 31, 32, 41, 42, 43, 44, 45, 46, \
47, 48, 49, 50, 51, 52, 53, 999, 255]
47, 48, 49, 50, 51, 52, 53, 999, 255, 128]
capc = [0, 0.18, 0.18, 0, 0.18, 0.18, 0, 1.69, 0, 0, 1.87, 3.03, 2.45, 0.18, 0, 0.18, 0.18, 0.18, 0.18, 0.18, 0, 0, \
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
return dict(zip(lceu, capc))
......