Commit c86b4e1a authored by roelo008's avatar roelo008

Testing met pollination model

Prep data for pest control
parent 85ba62d4
......@@ -16,7 +16,7 @@ VOORBEWERKINGEN
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
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
......@@ -43,6 +43,7 @@ bt_cats = [65535, 0, 101, 102, 103, 104, 201, 301, 401, 402, 403, 404, 501, 502,
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]
pest_control_cats_bt = []
bt_2_lceu = dict(zip(bt_cats, lceu_cats_bt))
# Dictionary Top10NL -> LCEU categorien
......
"""
Script for generating a LU map for the Pest Control model as developed by RIVM, March 2020. Map will have 5 LU cats:
PC-Cats
1: "eenjarig_bloemrijk"
2: "meerjarig_bloemrijk"
3: "meerjarig_bloemarm"
4: "bosjes"
5: "alleenstaande_bomenrij"
Map dimensions are equal to the LCEU map
Width x height: 27.000 x 32.000
Origin: 10.000, 620000 (RD-New)
Pixel size: 10,-1
Map is hierarchically composed of the following:
1. sloten_2019_10m_255.tif (sl)
2. Agrarisch Natuurbeheer (anb), source: S:\ESD_SpatialData\agrarisch_natuurbeheer\BiologischeLandbouw_RobSmidt 2019\ANLB123_union_NDB_TK
3. Beheertypenkaart 'huidig' (bt), source: BT_aangevuld_20190612time115357
4. Top10 (top10), source: TOP10NL_2018_seot_20181210time120347
Pre-processing
1. dit is extract van de Top10NL kaart, geresampled naar 10m via Aggregate, door Henk Meeuwsen. Ook reeds vertaald
naar PC-Cats
2: select relevante categorien uit field 'CODE_BEHEERPAKKET' en vertaal naar PC-Cats (result: ANB123_sel4pestcontrol.tif)
3: resample to 10m (NN) en export to TIF (result: BT_aangevuld_20190612time_resize.tif)
4. resample to 10m (NN) en export to TIF (result: TOP10NL_2018_sept_2018121_resize.tif)
In dit script
1. vertaal bt categorieren naar PC-Cats
2. vertaal top10 categorien naar PC-Cats
3. Voeg 1, 2, 3 en 4 hierarchisch samen in de volgorde:
Sloten \ ANB \ Beheertype \ Top10
Ie. gebruik 'Sloten' pixelwaardes, NoData opvullen met ANB pixelwaardes, resterende NoData opvullen met BT etc.
4. wegschrijven naar *.tif
Kaart gemaakt in opdracht Bart de Knegt, zie correspondentie 'Inputkaart voor plaagonderdrukking' maart 2020
Door: Hans Roelofsen, March 2020
"""
import os
import rasterio as rio
import numpy as np
from z_utils import esd
"""Lees bronbestanden als Rasterio object en als numpy array. Definieer ook NoData"""
sl = rio.open(r'c:\apps\temp_geodata\lceu_huidig2020\sloten_top10\Sloten_2019_10m_255.tif')
sl_arr = sl.read(1)
sl_nodata = 255
anb = rio.open(r'c:\apps\temp_geodata\lceu_huidig2020\ANB123_sel4pestcontrol.tif')
anb_arr = anb.read(1)
anb_nodata = 0
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
'''Map Beheertypen (bt) to Pest Control categorien. Note how NoData (65535) is mapped to 0'''
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]
pc_cats_bt = [0, 0, 0, 3, 3, 3, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 0, 0, 3, 0, 0, 3, 2, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2,
2, 0, 0, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 3, 3, 3,
0, 2, 0, 0, 0, 0, 0, 0, 2, 0, 2, 3, 3, 3, 2, 2, 2, 2, 2, 2, 3, 2, 4, 4, 3, 0, 3, 3, 3, 2, 3, 4, 4, 3, 4,
3, 3, 0, 0, 2, 4, 0, 0, 0, 0, 0, 0, 5, 0, 0, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 0, 0, 0, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 2, 2, 0, 4, 0, 0, 3,
4, 4, 4, 4, 4, 4, 4, 5, 4]
bt_2_pc = dict(zip(bt_cats, pc_cats_bt))
bt_pc_arr = esd.vec_translate(bt_arr, bt_2_pc)
'''Map Top10NL to Pest Control categorien. Note how NoData (32767) is mapped to 0'''
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]
pc_cats_top10 = [0, 3, 3, 0, 0, 0, 0, 0, 3, 0, 0, 3, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 0, 0,
3, 0, 0, 3, 0, 3, 3, 3, 0, 3, 0, 0, 3, 0, 3, 0, 0, 0, 0, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 0, 0, 0, 0,
0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 4, 0, 0, 0, 0, 3, 3, 3, 0, 0, 0, 0, 4, 4, 3, 0, 0, 0, 0, 0, 0, 0, 3, 3,
3]
top10_2_pc = dict(zip(top10_cats, pc_cats_top10))
top10_pc_arr = esd.vec_translate(top10_arr, top10_2_pc)
'''==OPTIE1=== Define where each array pops up in the output'''
# top10 where everyting else is NoData and Top10 self is not NoData
where_top10 = (bt_pc_arr == 0) & (anb_arr == 0) & (sl_arr == sl_nodata) & (top10_pc_arr != 0)
where_bt = (anb_arr == 0) & (sl_arr == sl_nodata) & (bt_pc_arr != 0) # bt where sloten AND anb is NoData
where_anb = (sl_arr == sl_nodata) & (anb_arr != 0) # anb where sloten is NoData
where_sl = sl_arr != sl_nodata # sloten where sloten is not NoData
out_arr = np.zeros_like(sl_arr).astype(np.uint8) # zeros array
out_arr[where_top10] = top10_pc_arr[where_top10] # put in top10 values where needed
out_arr[where_bt] = bt_pc_arr[where_bt] # put in bt values where needed
out_arr[where_anb] = anb_arr[where_anb] # put in anb values where needed
out_arr[where_sl] = sl_arr[where_sl] # put in sloten where needed
source_arr = np.zeros_like(sl_arr)
source_arr[where_sl] = 10 # alle cellen waar de waarde van sloten komt
source_arr[where_anb] = 20 # idem voor ANB
source_arr[where_bt] = 30 # idem voor BeheerType
source_arr[where_top10] = 40 # idem voor Top10
'''SAMPLE DATA
sloten = np.array(np.random.randint(low=0, high=4, size=49)).reshape(7, 7)
anb = np.array(np.random.randint(low=0, high=4, size=49)).reshape(7, 7)
bt = np.array(np.random.randint(low=0, high=4, size=49)).reshape(7, 7)
top10 = np.array(np.random.randint(low=0, high=4, size=49)).reshape(7, 7)
where_none = (bt == 0) & (sloten == 0) & (anb == 0) & (top10 == 0 )
where_top10 = (bt == 0) & (sloten == 0) & (anb == 0) & (top10 != 0)
where_bt = (anb == 0) & (sloten == 0) & (bt != 0)
where_anb = (sloten == 0) & (anb != 0)
where_sloten = (sloten != 0)
zeros = np.zeros_like(bt)
zeros[where_sloten] = 10
zeros[where_anb] = 20
zeros[where_bt] = 30
zeros[where_top10] = 40
zeros[where_none] = 0
'''
"""Write to file"""
out_dir = r'c:\apps\temp_geodata\lceu_huidig2020\out_pestcontrol_20200324'
prof = bt.profile
prof.update(dtype=rio.uint8)
with rio.open(os.path.join(out_dir, 'pest_control_land_use_20200324_v2.tif'), 'w', **prof) as dst:
dst.write(out_arr.astype(rio.uint8), 1)
with rio.open(os.path.join(out_dir, 'pest_control_land_use_20200324_sources.tif'),'w', **prof) as dst:
dst.write(source_arr.astype(rio.uint8), 1)
'''
Bestuivingsmodel
Theorie en idee: Marjolein Lof
Code: Hans Roelofsen
'''
import pandas as pd
import rasterio as rio
import geopandas as gp
import numpy as np
import scipy
from z_utils import esd
'''
INPUT VARIABLES
'''
# Basisregistratie Percelen
brp_src = r'c:\apps\temp_geodata\2G\brp15.tif' # RASTER
# LCEU
lceu_src = r'x:\ESD_SpatialData\EcosysteemUnits_LCEU\b_Huidig\versie_met _ANB_20191223\lceu_huidig_20191223.tif'
# Tabellen
hab_suitability_src = r'c:\apps\proj_code\esd\d_resources\2G\HabSuitPol.csv'
gewas_opbrengst_src = r'c:\apps\proj_code\esd\d_resources\2G\LUT_LEI_Eha_z.csv'
# Output
out_dir = r'c:\apps\proj_code\esd\c_out\2G'
'''
Read into memory
'''
brp = rio.open(brp_src) # rasterio object
lceu = rio.open(lceu_src) # rasterio object
brp_arr = brp.read(1) # numpy array
lceu_arr = lceu.read(1) # numpy array
hab_suit_tab = pd.read_csv(hab_suitability_src, sep=';', comment='#') # pandas dataframe
gew_opbr_tab = pd.read_csv(gewas_opbrengst_src, sep=';', comment='#') # pandas dataframe
'''
Check consistency in size
'''
esd.check_affine_consistency(brp, lceu)
'''
Calculate dependency of crops on pollination
'''
brp_afh_dict = dict(zip(gew_opbr_tab.brpcode, gew_opbr_tab.afhankelijkheid_bestuiving))
if set(np.unique(brp_arr).tolist()).difference(set(gew_opbr_tab.brpcode)):
print('hallo')
# todo: extend brp_afh_dict
bpr_afh = esd.vec_translate(brp_arr, brp_afh_dict) # perceel afhankelijkheid van bestuiving [5, 25, 64, 95]
'''
Calculate suitability of land cover types to shelter pollinators
'''
hab_suit_dict = dict(zip(hab_suit_tab.EUcode, hab_suit_tab.habsuit))
hab_suit = esd.vec_translate(lceu_arr, hab_suit_dict)
# spatial scale pollination NL
# See: C:\apps\proj_code\esd\b_analysis\2G_FFT.r
nr = 10000 # number of grid points in both x,y directions
dx, dy = 20, 20 # size of gridpoints m
dxkm, dykm = dx/1000, dy/1000 # size of gridpoints in km
xmin = -(np.divide(nr, 2))*dxkm # lowest X value
ymin = -(np.divide(nr, 2))*dxkm # lowest Y value
xlmin = xmin + dxkm/2 # minimum x coordinate, xmin + half a cell
ylmin = ymin + dykm/2 # minimum y coordinate, xmin + half a cell
# x- and y-coordinates [-99.99, -99.97, -99.95 ... 99.95, 99.97, 99.99]
x = np.arange(xlmin, np.multiply(np.divide(nr, 2), dxkm), dxkm)
y = np.arange(ylmin, np.multiply(np.divide(nr, 2), dykm), dykm)
Y, X = np.meshgrid(x, y) # two arrays with X and Y coordinates
kernel = np.exp(np.multiply(-0.53, np.sqrt(np.add(np.square(X), np.square(Y))))) # dispersal kernel
n_kernel = np.divide(kernel, np.multiply(np.product([dxkm, dykm]), np.sum(kernel))) # normalized kernel
f_n_kernel = np.fft.fft2(n_kernel) # fft kernel. Hier beginnen veranderingen tov R fft: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/fft.html
pollination_kernel = np.multiply(np.product([dxkm, dykm]), f_n_kernel)
{
"bt2lceu": {
"from": "BT_aangevuld_20190612time115357",
"to": "LCEU",
"map": {
}
}
}
\ No newline at end of file
......@@ -9,10 +9,10 @@ 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, \
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, 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]
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]
return dict(zip(lceu, capc))
......@@ -81,4 +81,4 @@ def check_affine_consistency(*rios):
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
raise Exception(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