Commit 6cdfa722 authored by Roelofsen, Hans's avatar Roelofsen, Hans

new script for translating stuff into LCEU categories and merging rasters

parent 1aa24ed1
......@@ -32,14 +32,16 @@ import os
import numpy as np
import rasterio as rio
from utils import esd
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\c_bau\LCEU10_bau_gdal.tif' # path to LCEU scenario
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
norm = 10 # norm voor concentratie
scen_naam = 'BAU' # scenario naam
px = 10 # pixel grootte in m
myear = esd.myear_tab() # invangcapaciteit per landgebruiks-categorie
nl_onshore = 33893 # area of onshore Netherlands in sq km
......@@ -72,9 +74,9 @@ true_perc = np.product([np.divide(true_area, nl_onshore), 100]) # divide by ref
'''Write outputs to file'''
# TODO
msg = 'For scenario {0} there are {1} pixels w value <= {2}, which equals to:' \
'{3} square km, which is {4}% of onshore area of the Netherlands ({5} km2)'\
.format(scen_naam, true_pxls, norm, true_area, true_perc, nl_onshore)
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:
f.write(msg)
......
"""
Helper script to generate a LCEU-flavour map hierarchically based on the following inputs
1 Agrarisch Natuurbeheer percelen:
~\ESD_SpatialData\agrarisch_natuurbeheer\BiologischeLandbouw_RobSmidt 2019\ANB_BioLB_2018.gdb\ANB123_union_ANLB_NDB_TK
2 Beheertypenkaart 'huidig' dd 20190612
~\ESD_SpatialData\Beheertypenkaart_MNP_NNN\a_HUIDIG\NVK_export.gdb\BT_aangevuld_20190612time115357
3 Top10NL 'huidig'
~\ESD_SpatialData\Top10-verrasterd_2.5m\BNL_export_topografie.gdb\TOP10NL_2018_sept_20181210time120347
VOORBEWERKINGEN
1: select relevante categorien uit field 'CODE_BEHEERPAKKET' (zie: https://sharepoint.wur.nl/sites/NVK/_layouts/15/WopiFrame.aspx?sourcedoc=/sites/NVK/Shared%20Documents/5.%20Data/vertaaltabel_LCEU_naar_BBNN_en_Top10_rondstuur.xlsx&action=default
, rasterize to 10m en set value naar '27'
2: resample to 10m (NN) en export to TIF
3. resample to 10m (NN) en export to TIF
In dit script:
vertaal 2 naar LCEU categorien, zie: https://sharepoint.wur.nl/sites/NVK/Shared%20Documents/5.%20Data/vertaaltabel_LCEU_naar_BBNN_en_Top10_rondstuur.xlsx?Web=1
vertaal 3 naar LCEU categorien, zie: https://sharepoint.wur.nl/sites/NVK/Shared%20Documents/5.%20Data/vertaaltabel_LCEU_naar_BBNN_en_Top10_rondstuur.xlsx?Web=1
voeg 1, 2, 3 samen
Dit bestand wordt gemaakt op verzoek BdK eind december 2019 tbv Ecosysteemdienst modellen.
"""
import os
import rasterio as rio
import numpy as np
from z_utils import esd
# Dictionary Beheertypen -> LCEU categorien
bt_cats = [65535, 0, 101, 102, 103, 104, 201, 301, 401, 402, 403, 404, 501, 502, 601, 602, 603, 604, 605, 606, 701, 702, 801, 802,
803, 804, 901, 1001, 1002, 1101, 1201, 1202, 1203, 1204, 1205, 1206, 1301, 1302, 1401, 1402, 1403, 1501,
1502, 1601, 1602, 1701, 1702, 1703, 1704, 15001, 15002, 15003, 15011, 15012, 15121, 15122, 15123, 15124,
15125, 15126, 15127, 15128, 15129, 15130, 15131, 15132, 15133, 15134, 15135, 15251, 15252, 15253, 15354,
15455, 1, 2, 10001, 10201, 10101, 10301, 15004, 80201, 80202, 80203, 80204, 80205, 80206, 80207, 80208,
80209, 80210, 80211, 80212, 80213, 80214, 50101, 50102, 50103, 50104, 50105, 50106, 50107, 50108, 50109,
50110, 50111, 50112, 50113, 50114, 80215, 50115, 15013, 15014, 15015, 0, 10002, 10003, 10004, 10005, 10006,
10007, 140101, 140102, 140103, 140201, 140202, 140203, 150101, 150102, 150103, 150201, 150202, 150203,
160101, 160102, 160103, 160201, 160202, 160203, 170301, 170302, 170303, 150104, 150105, 10202, 10203, 10204,
140106, 140107, 140206, 140204, 140205, 150204, 150205, 150206, 150207, 150208, 160104, 160105, 160106,
160107, 160108, 160205, 160206, 160208, 170305, 170306, 10102, 10401, 10402, 10403, 10501, 40201, 40202,
140109, 140102, 160207, 140208, 150108, 170308, 10601, 10701, 15501, 15502]
lceu_cats_bt = [0, 0, 51, 11, 31, 27, 53, 53, 52, 52, 51, 51, 26, 26, 26, 26, 24, 24, 52, 52, 24, 25, 13, 11, 11, 11, 32, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 4, 21, 21, 21, 22, 21, 22, 21, 21, 21, 21, 21, 4, 1, 4, 27, 27, 52, 21, 5, 21, 5, 5, 5, 21, 2, 5, 21, 21, 5, 26, 27, 28, 28, 28, 29, 29, 29, 29, 45, 4, 53, 53, 4, 11, 11, 11, 27, 11, 11, 27, 27, 24, 26, 24, 22, 21, 24, 26, 26, 26, 26, 26, 26, 21, 22, 26, 21, 26, 26, 25, 31, 27, 21, 4, 1, 1, 0, 29, 41, 5, 29, 1, 45, 21, 21, 23, 21, 22, 23, 21, 22, 23, 21, 22, 23, 21, 22, 23, 21, 22, 23, 21, 22, 23, 5, 21, 53, 53, 53, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 25, 27, 27, 1, 11, 52, 52, 27, 22, 21, 21, 11, 21, 2, 27, 2, 21]
bt_2_lceu = dict(zip(bt_cats, lceu_cats_bt))
# Dictionary Top10NL -> LCEU categorien
top10_cats = [32767, -14210, -14200, -14190, -14182, -14180, -14162, -14160, -14140, -14130, -14125, -14100, -14080, -14060, -14020, -14010, -14000, -13300, -13000, -12800, -12700, -12610, -12600, -12505, -12500, -12430, -12420, -12415, -12410, -12405, -12400, 10000, 10100, 10101, 10200, 10201, 10202, 10300, 10301, 10310, 10311, 10400, 10401, 10410, 10411, 10412, 10500, 10501, 10510, 10511, 10512, 10600, 10601, 10602, 10700, 10701, 10710, 10711, 10720, 10721, 10730, 10731, 10740, 10741, 10742, 10750, 10751, 10752, 10760, 10761, 10780, 10781, 12400, 12405, 12410, 12415, 12420, 12425, 12430, 12435, 12500, 12505, 12600, 12605, 12610, 12700, 12800, 12810, 12820, 13000, 13100, 13200, 13300, 13400, 14000, 14002, 14010, 14012, 14015, 14020, 14022, 14030, 14040, 14050, 14052, 14060, 14062, 14070, 14080, 14082, 14090, 14100, 14110, 14120, 14130, 14132, 14135, 14140, 14142, 14145, 14160, 14162, 14163, 14165, 14170, 14175, 14180, 14182, 14183, 14190, 14192, 14195, 14200, 14205, 14210, 14212, 14215]
lceu_cats_top10 = [0, 11, 29, 12, 45, 45, 29, 29, 24, 4, 2, 29, 21, 23, 45, 1, 45, 41, 41, 52, 51, 26, 26, 52, 52, 53, 53, 53, 53, 53, 53, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 53, 53, 53, 53, 53, 53, 53, 53, 52, 52, 26, 26, 26, 51, 52, 52, 52, 44, 44, 45, 41, 3, 45, 45, 1, 1, 1, 45, 45, 41, 2, 2, 2, 23, 23, 21, 21, 21, 22, 29, 29, 2, 4, 4, 4, 24, 24, 24, 29, 29, 29, 29, 21, 21, 45, 45, 45, 25, 25, 12, 29, 29, 11, 11, 11]
top10_2_lceu = dict(zip(top10_cats, lceu_cats_top10))
# Lees bronbestanden als Rasterio object en als numpy array. Definieer ook NoData
bt = rio.open(r'c:\apps\temp_geodata\lceu_huidig2020\BT_aangevuld_20190612time_resize.tif')
bt_arr = bt.read(1)
bt_nodata = 65535
top10 = rio.open(r'c:\apps\temp_geodata\lceu_huidig2020\TOP10NL_2018_sept_2018121_resize.tif')
top10_arr = top10.read(1)
top10_nodata = 32767
anb = rio.open(r'c:\apps\temp_geodata\lceu_huidig2020\ANB123_sel_resize.tif')
anb_arr = anb.read(1)
anb_nodata = 0
# Vertaal elk naar LCEU categorien (ANB is al in geldigde LCEU categorien).
bt_lceu = esd.vec_translate(bt_arr, bt_2_lceu)
top10_lceu = esd.vec_translate(top10_arr, top10_2_lceu)
# Merge BT data into ANB data where ANB <> 27. np.where(condition, x, y) where true, yield x, otherwise y
arr1 = np.where(anb_arr == 0, bt_lceu, anb_arr)
# Merge Top10 data into remaining NoData (0) pixels of arr1
arr2 = np.where(arr1 == 0, top10_lceu, arr1)
# Write to file.
prof = bt.profile
with rio.open(r'c:\apps\temp_geodata\lceu_huidig2020\out\lceu_huidig_20191223.tif', 'w', **prof) as dst:
dst.write(arr2.astype(rio.uint16), 1)
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