Opbrengstenderving

Opbrengstdervingskaarten (WWL) in RuimteScanner

GeoDMS-locatie: /SourceData/Landbouw/Yield_Reduction GitHub-issue: ObjectVision/RSopen #500

Doel

Landsdekkende kaarten van de relatieve opbrengstderving in de landbouw als gevolg van de grondwatersituatie, op basis van de Waterwijzer Landbouw (WWL). Per cel een fractie tussen 0 en 1, afgeleid uit de huidige bodem, het huidige gewas, een gekozen beregeningsscenario en de grondwaterstand (GxG) van een gekozen variant en zichtjaar.

Invoer

De WWL-relatiedatabase, geleverd door Deltares als WWL_relatie_database.nc. Locatie: %RS_Lb_DataDir%/Gewassen/WWL_relatie_database.nc.

De bodemkaart BOFEK2020 op 25m (directe vergridding van het shapefile). Deze codering is gelijk aan de WWL-bodemcode, dus er is geen aparte vertaaltabel nodig.

De gewaskaart LGN op 25m (modus uit het 5m-bronbestand), via een reclass LGN naar WWL-gewascode (Classifications/Grondgebruik/LGNKlasse/WWL_ID), obv WWL documentatie.

De grondwaterstandskaarten GHG en GLG uit de Basisprognoses 2018 (beschreven in de Geactualiseerde Knelpuntenanalyse, Mens et al. 2019), per variant en zichtjaar, op 250m resolutie. 2050 gebruikt de meteoreeks 1929-2011, 2085 gebruikt de kortere reeks 1974-2003. Die laatste maakt 2050 en 2085 niet helemaal 1-op-1 vergelijkbaar; dat is een bewuste eigenschap van de bronkaarten.

De WWL-relatiedatabase verklaard

De nc is geen geografische kaart en geen platte tabel, maar een responssurface. Hij bevat vijf schadevariabelen: dmgwet (natschade), dmgdry (droogteschade), dmgdir (directe schade), dmgind (indirecte schade) en dmgtot (totale schade). Voor de hoofdindicator gebruiken we dmgtot. De waarden zijn schadepercentages (0 tot ongeveer 99), waarbij laag gezond betekent.

Elke variabele heeft 2844 banden = 79 bodems maal 9 gewassen maal 2 beregeningsopties maal 2 klimaatopties. Elke band is een grid van 301 bij 301. De twee rasterassen zijn GHG en GLG, niet de ruimte. Per cel staat het schadepercentage voor die combinatie van bodem, gewas, beregening en klimaat, bij die specifieke GHG en GLG.

De assen lopen van 0 tot 300 cm beneden maaiveld, met de relatie cm = 300 - pixelindex. De rij is de GHG-as, de kolom is de GLG-as. Alleen de helft waar GLG dieper is dan GHG is geldig; de andere helft is NaN, want de laagste grondwaterstand kan niet ondieper zijn dan de hoogste.

De as-toewijzing en de richting zijn empirisch vastgesteld: de natste hoek (maximale natschade) ligt op pixelindex (300, 300) en hoort bij GHG en GLG gelijk 0; de droogste hoek (maximale droogteschade) ligt op (0, 0) en hoort bij GHG en GLG gelijk 300. Een controlecel buiten de diagonaal bevestigde dat de rij de GHG-as is en de kolom de GLG-as.

Codes en dimensies

Bodem: 79 codes, 1001 tot en met 5007, oplopend. Dit is exact de nc-volgorde, gelijk aan de oplopende BOFEK_K-classificatie.

Gewas: 9 codes, namelijk 1, 6, 7, 9, 12, 13, 20, 22, 23. Let op: de volledige WWLKlasse heeft 23 gewassen; alleen deze 9 zitten in de nc en alleen die mogen in de lijst staan, in deze volgorde.

Beregening (irrigatie_optie): 2 opties. 0 is geen beregening, 1 is wel beregening. Geen gradatie in hoeveelheid.

Klimaat (zichtjaar_klimaat): 2 opties. 0 is de periode 1991-2020 (referentieklimaat), 1 is de periode 2036-2065 in het warme KNMI-scenario W-hoog. De toekomstige zichtjaren 2050 en 2085 worden gekoppeld aan klimaat 1, omdat de GxG-kaart de grondwaterkant levert en de WWL-klimaat de atmosferische forcering; een toekomstige GxG met referentieklimaat zou inconsistent zijn.

Rekenmethode

De banden worden eenmalig ingelezen met for_each_ndvs, per schadevariabele 2844 losse float32-attributen (b1 tot en met b2844), elk via BANDS=.

De bandindex wordt ontleed naar bodem, gewas, beregening en klimaat. De volgorde van binnen naar buiten is klimaat, beregening, gewas, bodem:

band      = id + 1
klimaat   =  id        % 2
irrigatie = (id /  2)  % 2
gewas     = (id /  4)  % 9
bodem     =  id / 36

De positie 0 tot en met 8 voor gewas en 0 tot en met 78 voor bodem indexeert in de geordende codelijsten hierboven.

Per variant, zichtjaar en beregening gebeurt dan het volgende, in de template OpbrengstenDerving_T:

Selecteer met select_with_org_rel de 711 banden die bij de gekozen beregening en klimaat horen (combo_sel).

Bouw lut = combine(combo_sel, wwl_grid), een tabel van 711 maal 90601. Per rij wordt de GHG en GLG in cm afgeleid uit de pixelpositie (ghg = 300 - rij, glg = 300 - kolom) en wordt een sleutel gemaakt uit bodem, gewas, ghg en glg.

Stapel de 711 banden in combo_sel-volgorde tot een enkel attribuut over lut met union_data. De argumentenlijst wordt gegenereerd met AsItemList en uitgevoerd via een indirecte expressie. Doordat lut combo_sel als buitenste as heeft, valt de gestapelde array precies op de juiste rijen:

attribute<float32> dmgtot_stack := ='union_data(., ' + AsItemList('WWL_relatie_database/dmgtot_bands/b' + string(combo_sel/band)) + ')';

De sleutel codeert de vier velden collisievrij in een uint64 (32 is groter dan de hoogste gewascode, 512 is groter dan 300):

key = ((bodem_code * 32 + gewas_code) * 512 + ghg) * 512 + glg

Landsdekkend (op AdminDomain, 25m) wordt per cel de bodemcode (BOFEK), de gewascode (LGN naar WWL) en de GxG-index bepaald. De GxG-index is de grondwaterstand in cm, truncated op 0 tot 300, met nulls behouden:

ghg_idx = IsDefined(ghg_cm) ? int32(max_elem(0f, min_elem(300f, ghg_cm)) + 0.5f) : null_i

Met dezelfde sleutel wordt per cel een rlookup op lut/key gedaan, en de gevonden dmgtot_stack gedeeld door 100 levert de opbrengstderving als fractie.

Parameters en aannames

gxg_sign = 1: positief in de GxG-kaart betekent diepte onder maaiveld. Geverifieerd met de IJsselmeer-cellen (negatief, water boven maaiveld) en de Veluwe (positief, diep).

gxg_scale = 100: de GxG-kaarten staan in meters en worden naar cm omgezet. Bevestigd doordat de derving alleen met deze schaal fysiek plausibel werd.

zichtjaar_klimaat = 1 voor alle toekomstige zichtjaren.

dmgtot is schade, niet relatieve opbrengst. Bevestigd doordat het landelijk gemiddelde rond 0,10 ligt met het zwaartepunt onderaan; bij relatieve opbrengst zou dat rond 0,9 liggen.

Varianten en uitvoer

De varianten worden vermenigvuldig NL2120_Varianten (BAU1, BAU2), WWL_Zichtjaren (2050, 2085) en WWL_BeregeningsK (GeenBeregening, WelBeregening). Dat geeft 8 kaarten, allemaal bij klimaat 1.

De template OpbrengstenDerving_T wordt via for_each_ne voor elke combinatie geinstantieerd (container OpbrengstenDerving_Dynamisch). De GxG-paden worden dynamisch opgebouwd uit de variant- en zichtjaarnaam.

De uitvoer wordt per combinatie weggeschreven als GeoTIFF:

%RS_Lb_DataDir%/Gewassen/WWL_Opbrengstderving/<variant>_<zichtjaar>_<beregening>.tif

Bekende beperkingen

De GxG staat op 250m, dus de grondwaterkant van de derving is blokkerig op 250m, terwijl bodem en gewas op 25m varieren.

Ongeveer 13 miljoen landbouwcellen (rond 8000 km2) hebben geen GxG-dekking en krijgen dus geen derving. Dit dekkingsgat moet nog onderzocht en verantwoord worden.

2085 gebruikt een kortere meteoreeks dan 2050, dus die twee zijn niet helemaal 1-op-1 vergelijkbaar.

De componenten natschade, droogteschade, directe en indirecte schade tellen niet noodzakelijk op tot de totale schade; gebruik dmgtot als hoofdindicator.