Opbrengstenderving

GeoDMS-locatie: /SourceData/Landbouw/Yield_Reduction

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 pagina beschrijft de technische uitwerking. Wat droogteschade en natschade zijn, uit welke bron de relaties komen en welke beperkingen van de Waterwijzer Landbouw daarbij gelden, staat op Droogteschade en natschade.

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.

Het beregeningsregime per cel. Beregening is geen aparte kaartdimensie meer maar zit per cel in de sleutel, en welke waarde een cel krijgt volgt uit ModelParameters/Landbouw/BeregeningScenario_name: A geeft nergens beregening, B1 overal, en B2 volgt de kaart Beregeningskaart_Basisprognoses_2018 op AdminDomain, dus alleen de percelen die in 2018 een installatie hadden. De default is B2, en dat is de afspraak uit #542.

Tot #780 werd die parameter niet gelezen: de actieve route las onvoorwaardelijk de kaart van 2018, terwijl de default op B1 stond. Het model rekende dus B2 en de parameter zei B1. Wie de uitkomsten van voor die koppeling wil reproduceren zet het scenario op B2. Van B2 naar B1 verlaagt de derving op landbouwgrond met 0,6 tot 5,5 procentpunt, gemeten op zichtjaar 2040 van BAU.

De grondwaterstandskaarten GHG en GLG uit de Basisprognoses 2018 (beschreven in de Geactualiseerde Knelpuntenanalyse, Mens et al. 2019), per hydrologische levering en zichtjaar, op 250m resolutie. Er zijn drie leveringen: BAU, BAU2 en, sinds 8 september 2026, NbSGenuanceerd. Die derde is geen aparte som maar de BAU-stand met daaroverheen de grondwaterstand van de veenbouwsteen die op de cel is gealloceerd (#263, #782). Welke variant welke levering leest staat in VariantParameters/VariantK/Hydrologie_Levering; meer dan een variant mag dezelfde levering aanwijzen. 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 (beregening in WWL_relatie_database/combo): 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.

Welke van de twee een kaart krijgt is sinds #782 een scenariokeuze en geen variantkeuze, en staat in ModelParameters/Landbouw/Scenario. Die container heeft twee helften: Klimaat_name, die op Warm staat, en Groei_name, die volgt uit het scenario dat de casus draait. Samen wijzen zij de band WarmHoog aan. De kaart van 2017 leest de referentieband, want daar staan waargenomen grondwaterstanden tegenover; de kaarten van 2050 en 2085 lezen de scenarioband.

De twee helften zijn niet los in te stellen. De database kent alleen de banden Referentie en WarmHoog, dus een combinatie als WarmLaag bestaat er niet en zou op een lege verwijzing uitkomen. De groeihelft draagt in deze keten dus geen eigen informatie; zij zit alleen in de naam van de band.

Tot 8 september 2026 stond hier een tak op de variantnaam, die BAU op WarmHoog zette en BAU2 op Referentie. Deltares heeft in #782 vastgelegd dat het klimaat geen variantkenmerk is en dat deze studie in alle varianten naar Warm kijkt. Daarmee komt wel een tegenspraak terug: BAU2 leest de referentiestanden R2050BP18 en R2085BP18 en krijgt nu weer de warme schaderelatie. Dat is een bewuste keuze van Deltares en geen omissie, maar of BAU2 daarmee nog de bedoelde gevoeligheidsvariant is staat als vraag open in #782.

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 GDAL-optie BANDS is eengebaseerd: het geldige bereik is 1 tot en met 2844, en b<N> moet dus band N vragen.

Daar zat tot 8 september 2026 een fout in (#263). De optie werd gevuld met de nulgebaseerde rij-index, dus b1 vroeg band 0 en viel om, en elk ander item bN las band N-1, oftewel de band van de combinatie ernaast. Omdat klimaat de snelst wisselende as is, kreeg elke WarmHoog-aanvraag de referentiewaarde van dezelfde bodem, hetzelfde gewas en dezelfde beregening terug, en een referentieaanvraag de WarmHoog-waarde van een beregeningsstap lager. Gemeten op gras zonder beregening op klei: 18,09 in plaats van 18,46, en met beregening 18,46 in plaats van 16,61. De fout kwam pas boven toen om de basisjaarkaart van gras werd gevraagd, want die vraagt als enige band 1 en viel daarmee hard om; alle andere combinaties gaven een geldig ogend getal van de verkeerde rij, en het proces liep af met exitcode 0.

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
beregening = (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 volgt ModelParameters/Landbouw/Scenario: de referentieband voor de kaart van 2017, de scenarioband WarmHoog voor 2050 en 2085.

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, NbSGenuanceerd), WWL_Zichtjaren (2017, 2050, 2085) en de 9 WWL_Gewassen. Dat geeft 3 maal 3 maal 9 is 81 kaarten. NL2120_Varianten is ondanks zijn naam geen variantenlijst maar de lijst met hydrologische leveringen; welke variant welke leest staat in VariantK. 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:

%LocalDataProjDir%/BaseData/Landbouw/WWL_Opbrengstderving/<levering>_<zichtjaar>_<gewas>_<studiegebied>.tif

Dat pad stond tot 5 september 2026 onder %RS_Lb_DataDir%, dus in de gedeelde brondata, en werd bij elke run herschreven. Dit is geen brondata: de invoer is de relatiedatabase plus de grondwaterstanden, en die staan er los naast. De suffix met het studiegebied hoort erbij omdat AdminDomain met het studiegebied meebeweegt; zonder die suffix leest een deelgebiedrun het landsdekkende bestand in een kleiner domein, en dat geeft een lege uitkomst zonder foutmelding.

Het maken van de kaarten is de zware kant van de keten, want per gewas worden honderden banden uit de netcdf gelezen. Daarom leest de rest van het model standaard de weggeschreven tif terug en niet het rekenende item; de schakelaar is ModelParameters/OpbrengstdervingOntkoppeld. Wie de kaarten opnieuw wil maken vraagt de haak OpbrengstenDerving_Dynamisch/generate op, die met opzet de schrijfkant aanroept en niet de leeskant.

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.