Opbrengstenderving

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 gewas (waarvoor de kaart wordt gemaakt), het beregeningsregime per cel en de grondwaterstand (GxG) van een gekozen variant en zichtjaar.

Deze derving wordt dynamisch per variant, zichtjaar en gewas uitgerekend en gaat als opbrengstderving rechtstreeks de landbouw-suitability in (zie Landbouw methode, sectie Opbrengstderving). Dit verving de eerdere, statisch ingelezen WWL2-dervingskaarten per beregeningsscenario (issue #539). De opzet is uitgewerkt in issue #500.

Invoer

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

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 wordt niet als kaart ingelezen: per dervingskaart is het gewas een vaste parameter (Gewas_code). Welke LGN-klasse naar welke WWL-gewascode gaat, is vastgelegd in Classifications/Grondgebruik/LGNKlasse/WWL_ID, obv WWL documentatie; die reclass wordt in de landbouwmodule gebruikt om per cel het juiste gewas op te zoeken.

De beregeningskaart Beregeningskaart_Basisprognoses_2018 op AdminDomain. Hieruit volgt per cel of er wel of geen beregening is (WWL_Beregenings_rel). Beregening is dus geen aparte kaartdimensie meer, maar wordt per cel uit deze kaart bepaald.

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. In de huidige implementatie wordt voor alle dervingskaarten klimaat 1 (W-hoog) gebruikt (Zichtjaar_klimaat := WWL_Klimaatscenarios/v/WarmHoog), omdat de GxG-kaart de grondwaterkant levert en het WWL-klimaat de atmosferische forcering; de combinatie met een toekomstige GxG en referentieklimaat zou inconsistent zijn. De klimaatkeuze is daarmee losgekoppeld van het zichtjaar.

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 gewas gebeurt dan het volgende, in de template OpbrengstenDerving_T:

Selecteer met select_with_org_rel de 158 banden die bij het gekozen gewas en klimaat horen (combo_sel). Dat zijn 79 bodems maal 2 beregeningsopties; beide beregeningsopties blijven dus in de tabel, want de keuze daartussen valt pas per cel. (De oude per-beregening-opzet selecteerde 711 banden per beregening; dat is vervangen.)

Bouw lut = combine(combo_sel, wwl_grid), een tabel van 158 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 beregening, bodem, gewas, ghg en glg.

Stapel de 158 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 vijf velden collisievrij in een uint64 als mixed-radix getal (slots: beregening < 2, bodem_code < 8192 omdat de hoogste BOFEK_K-code 5007 is, gewas < 32, ghg < 512, glg < 512):

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

Landsdekkend (op AdminDomain, 25m) wordt per cel de bodemcode (BOFEK), de beregeningscode (uit de beregeningskaart), de gewascode (vast per kaart) 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

Per cel wordt met exact dezelfde sleutelformule een rlookup op lut/key gedaan; doordat de beregeningscode in de sleutel zit, selecteert elke cel automatisch de juiste van de twee beregeningsbanden. 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 uitgevraagd met combine_uint8 over drie assen: NL2120_Varianten (BAU, BAU2), WWL_Zichtjaren (2017, 2050, 2085) en de 9 WWL_Gewassen. Dat geeft 2 maal 3 maal 9 = 54 kaarten, allemaal bij klimaat 1 (W-hoog). Beregening is geen as meer (die zit per cel in de sleutel), in tegenstelling tot de oude opzet met WWL_BeregeningsK als dimensie.

De template OpbrengstenDerving_T wordt via for_each_ne voor elke combinatie geinstantieerd (container OpbrengstenDerving_Dynamisch, met onderliggende padstructuur <variant>/<zichtjaar>/<gewas>). 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>_<gewas>.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.