Commit 53839a54 authored by roelo008's avatar roelo008

bevolkingsraster toegevoegd en rapportage over # inwoners wonende onder en boven de norm

parent 22c019b3
......@@ -13,12 +13,14 @@ LCEU landgebruikskaart volgens scenario, 10m resolutie (lceu_scen, raster)
MYEAR: afvangcapaciteit per landgebruikstype
norm-waarde concentratie (norm, float)
scenario naam (naam_scen, string)
CBS 100m vierkant shapefile met bewonersaantallen per hok]\
Output
1. Kaart van concentratie zoals die zou zijn in afwezigheid van vegetatie (conc_blanco)
2. Kaart van concentratie zoals die verwacht wordt in het scenario (conc_scen)
3. Kaart als 2 boven/onder norm (eval_scen)
4. areaal onder/boven de norm (score_scenario)
5. aantal inwoners onder de norm/aantal inwoners boven de norm.
Berekeniningen
conc_blanco = conc_rivm * (1 + (myear * 0.137)/(1-myear*0.137)) (note output >= input)
......@@ -30,6 +32,7 @@ input rasters zijn geharmoniseerd, ie identieke Affine
'''
import os
import datetime
import numpy as np
import rasterio as rio
......@@ -38,22 +41,30 @@ from z_utils import esd
'''INPUT'''
# 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\e_hdbb\LCEU10_hdbb_gdal.tif' # path to LCEU scenario
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
scen_naam = 'hdbb' # scenario naam
norm = 10 # norm voor concentratie
px = 10 # pixel grootte in m
myear = esd.myear_tab() # invangcapaciteit per landgebruiks-categorie
px = 10 # pixel grootte in m van de LCEU kaart en eigenlijk ook van alle andere rasters.
nl_onshore = 33893 # area of onshore Netherlands in sq km
myear = esd.myear_tab() # invangcapaciteit per landgebruiks-categorie
'''Read rasters into memory'''
conc_rivm = rio.open(conc_rivm_in) # rasterio raster object
lceu_orig = rio.open(lceu_orig_in)
lceu_scen = rio.open(lceu_scen_in)
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
'''Calculate blanco concentratie '''
myear_blanco = esd.vec_translate(a=lceu_orig_arr, my_dict=myear)
......@@ -66,31 +77,62 @@ conc_scen = esd.conc_scenario(conc_blanco, myear_scen)
'''Treshold conc_scen to the norm'''
eval_scen = np.where(conc_scen <= norm, 1, 0).astype(np.uint8) # set to 1 where value <= norm
'''Calculate population living below/above norm'''
pop_below_th = np.where(eval_scen == 1, cbs_bevol_arr, 0)
pop_above_th = np.where(eval_scen == 0, cbs_bevol_arr, 0)
'''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]) # divide by reference area
total_pop = cbs_bevol_arr.sum().astype(np.uint64)
below_pop = pop_below_th.sum().astype(np.uint64)
above_pop = pop_above_th.sum().astype(np.uint64)
'''Write outputs to file'''
prof = lceu_orig.profile
out_dir = r'.\c_out'
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'
out_dir = os.path.join(out_base_dir, '{0}{1}'.format(timestamp_brief, scen_naam))
os.mkdir(out_dir)
basename = 'ESD2J_rpm25_{0}'.format(scen_naam)
prof = lceu_orig.profile
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
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
prof.update(dtype=rio.float32)
with rio.open(os.path.join(out_dir, basename+'_pop_below_tresh.tif'), 'w', **prof) as dst:
dst.write(pop_below_th, 1) # write populating living below treshold raster
with rio.open(os.path.join(out_dir, basename+'_pop_above_tresh.tif'), 'w', **prof) as dst:
dst.write(pop_above_th, 1) # write population living above treshold raster
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)
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
'''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'.\c_out\{0}_report.txt'.format(os.path.splitext(scen_naam)[0]), 'w') as f:
f.write(msg)
footer = '\nMade with Python 3.5 using pandas, geopandas, by Hans Roelofsen, WEnR team B&B,' \
' dated {0}'.format(timestamp_full)
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, os.path.basename(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} live in areas where ' \
'the treshold is exceeded and {2} in areas below the treshold.'.format(total_pop, above_pop, below_pop)
print(msg1)
print(msg2)
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(footer)
......
'''Helper functions for several ESD models'''
import os
import rasterio as rio
import numpy as np
def myear_tab():
......@@ -53,3 +55,26 @@ def vec_translate(a, my_dict):
https://stackoverflow.com/questions/16992713/translate-every-element-in-numpy-array-according-to-key
'''
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'''
# unpack arguments
rio_lst = [x for x in rios]
# list of affine transformations of each
try:
affs = [getattr(x, 'affine') for x in rio_lst]
except AttributeError:
raise Exception('A non-rasterio object was presented for Affine inspection')
# check for equality of all elements in the list
if affs.count(affs[0]) == len(affs):
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
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