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) van WUR. Per cel een fractie tussen 0 en 1, afgeleid uit de bodem, het gewas waarvoor de kaart wordt gemaakt, het beregeningsregime per cel en de grondwaterstand (GHG en GLG) van een hydrologische levering in een ankerjaar.

De kaarten gaan als opbrengstderving rechtstreeks de landbouwgeschiktheid in. Hoe een zichtjaar zijn derving uit de ankerjaren afleidt staat onderaan de rekenmethode; hoe de derving in het saldo per hectare doorwerkt staat op Landbouw methode, paragraaf Opbrengstderving.

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.

Invoer

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

De bodemkaart BOFEK2020, als shapefile direct vergrid naar AdminDomain (SourceData/Bodem/BOFEK/BOFEK_CODE). De BOFEK-code is dezelfde als de bodemcode van de WWL, dus er is geen vertaaltabel nodig. De shapefile heeft 5.771 polygonen zonder code, samen een miljoen hectare water, bebouwing en buitenland, met daarin 35.000 ha landbouw van het basisjaar: kassen die BOFEK als bebouwd ziet, oeverstroken en dorpsranden. Die cellen, en de cellen buiten de kaart, krijgen de code van de dichtstbijzijnde gecodeerde cel van 250 meter (BOFEK_CODE_Gevuld); een cel van 25 meter zonder code in een gecodeerde cel van 250 meter krijgt die code direct. Welke cellen gevuld zijn staat in BOFEK_CODE_IsAangevuld.

De gewaskaart wordt niet ingelezen: per dervingskaart is het gewas een vaste parameter (Gewas_code); de landbouwmodule rekent per gewas een geschiktheid en zoekt niet per cel een gewas op. De koppeling met de gewassen van het model loopt via de kolom GeneriekeNaam in Classifications/Grondgebruik/WWLKlasse, die letterlijk gelijk moet zijn aan de gewasnaam in Classifications/Landbouw/GewasSoortYR: Templates/Landbouw.dms zoekt de dervingskaart onder die naam op, en een WWL-gewas zonder generieke naam krijgt wel een kaart maar niets leest die. De tabel van LGN-klasse naar WWL-gewascode in Classifications/Grondgebruik/LGNKlasse/WWL_ID (WWL-documentatie) leest de dervingsketen niet.

Het beregeningsregime per cel. Beregening is geen as van de kaartenset maar zit per cel in de opzoeksleutel, 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 van Deltares, dus alleen de percelen die in 2018 een installatie hadden (de codes 1 en 2 van die kaart). B2 is de standaardstand, de afspraak uit #542; B1 verlaagt de derving op landbouwgrond met enkele procentpunten. Op de cellen waar een hydrologische levering de grondwaterstand van een veenbouwsteen oplegt staat de beregening uit, wat het scenario ook zegt, want bij natte teelt en bij de veenvormende bouwstenen hoort geen beregening (#542, #793). De regel raakt de WWL-gewassen op die cellen, en daarmee of de allocatie daar iets anders dan natte teelt zou willen neerzetten; riet en lisdodde lezen geen WWL-kaart.

De grondwaterstandskaarten GHG en GLG uit de Basisprognoses 2018 (beschreven in de Geactualiseerde Knelpuntenanalyse van Deltares, 2019), per hydrologische levering en ankerjaar, op 250 m, in SourceData/Water/Grondwaterstanden. Elke tif wordt op zijn eigen grid van 250 meter gelezen en de cellen zonder waarde krijgen de waarde van de dichtstbijzijnde gevulde cel (Grondwaterstanden/Kaart_T, met IsAangevuld als masker); de kale levering blijft beschikbaar als GLG_2017_TR29_11_Bron voor lezers die null als ontbrekende kennis willen zien, zoals de veenreserve in de koolstofboekhouding. Een levering is een eigen hydrologische som van Deltares, of een bestaande stand met daaroverheen de grondwaterstand die de gealloceerde veenbouwsteen per cel oplegt (Opgelegd), zodat de vernatting van een variant terugkomt in de derving van precies die cellen; een variant die de standen van een andere levering leest ziet haar eigen vernatting niet. Waar de bron van een bouwsteen alleen een GLG geeft schuift de GHG met de GLG mee, geklemd op maaiveld; of dat de bedoelde toestand is staat als vraag in #782. Welke variant welke levering leest staat in VariantParameters/VariantK/Hydrologie_Levering; meer varianten mogen dezelfde levering aanwijzen, en welke leveringen er zijn staat op Toepassing NL2120.

Het ankerjaar 2017 draagt voor elke levering de waargenomen referentiestand (REF2017BP18), en bij een levering met veenoplegging ook de opgelegde standen, want de drie kaarten zijn klimaatankers en geen beleidsjaren en de vernatting geldt vanaf het eerste zichtjaar volledig. 2050 gebruikt de meteoreeks 1929-2011, 2085 de kortere reeks 1974-2003; dat maakt de twee niet helemaal een op een vergelijkbaar en is een 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). Het model gebruikt 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. In de netcdf is de rij de GHG-as en de kolom de GLG-as, beide oplopend van 0 naar 300. GDAL leest de rij-as gespiegeld en de kolom-as niet, dus in de grid zoals het model hem ziet geldt ghg = 300 - rij en glg = kolom. 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 database documenteert de astoewijzing en de richting niet zelf; zij volgen uit de ligging van de twee uiterste hoeken in de grid zoals het model hem leest. De natste hoek, 99 procent schade bij GHG en GLG gelijk 0, ligt op rij 300 kolom 0; de droogste hoek, GHG en GLG gelijk 300, ligt op rij 0 kolom 300; en de onmogelijke helft is NaN waar de kolom kleiner is dan 300 min de rij. Die drie samen leggen beide assen en beide richtingen vast (#837).

Codes en dimensies

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

Gewas: 9 codes, namelijk 1, 6, 7, 9, 12, 13, 20, 22, 23 (WWL_relatie_database/WWL_Gewassen). 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 (klimaat in combo, per kaart 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 Classifications/Landbouw/WWL_Klimaatscenarios heten zij Referentie en WarmHoog.

Welke van de twee een kaart krijgt is een scenariokeuze en geen variantkeuze, en staat in ModelParameters/Landbouw/Scenario. Die container heeft twee helften: Klimaat_name, een instelling, en Groei_name, die volgt uit het scenario dat de casus draait. Samen vormen zij de naam van de band (WarmHoog), en een IntegrityCheck eist dat die naam in de database bestaat. Het ankerjaar 2017 leest altijd de referentieband (Basisjaar_name), want daar staan waargenomen grondwaterstanden tegenover; de kaarten van 2050 en 2085 lezen de scenarioband. Welke waarde de klimaathelft heeft is een keuze van de toepassing; in NL2120 staat hij in alle varianten op Warm, een afspraak met Deltares.

De twee helften zijn niet los in te stellen: de database kent alleen Referentie en WarmHoog, dus een combinatie als WarmLaag valt op de IntegrityCheck om, en de groeihelft draagt in deze keten geen eigen informatie.

Doordat het klimaat van de schaderelatie aan het scenario hangt en de grondwaterstanden aan de levering van de variant, krijgt een variant die referentiestanden leest toch de warme schaderelatie. Dat is een bewuste keuze van Deltares, maar of zo’n variant daarmee nog de bedoelde gevoeligheidsvariant is staat als vraag open in #782, net als de vraag of de standen van 2050 met de warme schaderelatie voor een zichtjaar als 2040 de bedoeling zijn.

Rekenmethode

De banden van dmgtot worden eenmalig ingelezen met for_each_ndvs, als 2844 losse float32-attributen (b1 tot en met b2844), elk via de GDAL-optie BANDS=. Die optie is eengebaseerd: het geldige bereik is 1 tot en met 2844, en b<N> vraagt band N, dus de rij-index van combo plus een. De containers dmgwet_bands en dmgdry_bands hebben dezelfde opbouw en worden door niets gelezen.

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 levering, ankerjaar 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 in de tabel, want de keuze daartussen valt pas per cel.

Bouw LUT = combine(combo_sel, wwl_k), een tabel van 158 maal 90601, waarin wwl_k de grid als platte reeks in rijvolgorde is (element k is rij k div 301 en kolom k mod 301). Per rij wordt de GHG en GLG in cm afgeleid uit de positie (ghg = 300 - rij, glg = 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_opslag = combine(combo_sel, dmgtot) met union_data. De argumentenlijst wordt gegenereerd met AsItemList en uitgevoerd via een indirecte expressie. Die stapel staat per band in de opslagvolgorde van de grid, en dat is de tegelvolgorde en niet de rijvolgorde; wwl_k/perm zegt per punt waar zijn waarde in die volgorde staat, en LUT/opslag_rel leest de stapel daarmee op de juiste positie. Zonder die permutatie hoort de gestapelde waarde bij een ander punt dan de sleutel zegt:

attribute<float32>    dmgtot_stack (LUT_opslag) := ='union_data(., ' + AsItemList('WWL_relatie_database/dmgtot_bands/b' + string(combo_sel/band)) + ')';
attribute<LUT_opslag> opslag_rel   (LUT)        := value(uint32(first_rel) * 90601 + uint32(wwl_k/perm[second_rel]), LUT_opslag);
attribute<float32>    dmgtot_stack (LUT)        := LUT_opslag/dmgtot_stack[opslag_rel];

De sleutel codeert de vijf velden collisievrij in een uint64 als mixed-radix getal (slots: beregening < 2, bodem_code < 8192 omdat de hoogste BOFEK-code 5007 is, gewas < 32, ghg < 512, glg < 512):

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

Landsdekkend, op AdminDomain, wordt per cel de bodemcode (BOFEK), de beregeningscode (scenario of kaart, met de uitzondering voor opgelegde vernatting), de gewascode (vast per kaart) en de GxG-index bepaald. De GxG-index is de grondwaterstand in cm, geklemd op 0 tot 300 en afgerond op hele centimeters, met nulls behouden:

ghg_index = 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. Door de vulling van de bodemcode en de grondwaterstanden heeft elke cel een sleutel; alleen een combinatie die in de relatiedatabase zelf ontbreekt zou nog een lege waarde geven, en de database is op de fysisch mogelijke helft volledig gevuld.

Van ankerjaar naar zichtjaar

De kaarten bestaan alleen voor de drie ankerjaren 2017, 2050 en 2085 (Classifications/Landbouw/WWL_Zichtjaren). De landbouwgeschiktheid van een zichtjaar (Templates/Landbouw/Suitability_T) leest per gewas de kaarten van de twee ankers die het zichtjaar omsluiten (WWL_Anker_Onder, WWL_Anker_Boven) en weegt ze lineair naar de afstand tot elk anker; WWL_Gewicht is het gewicht van het bovenste anker, geklemd tussen 0 en 1, zodat een zichtjaar op een anker alleen dat anker leest en een zichtjaar voorbij 2085 alleen de kaarten van 2085. Zo verloopt de derving geleidelijk over de reeks zichtjaren in plaats van in een sprong. Waar maar een anker een waarde heeft telt die ene (WWL_dyn); waar geen van beide iets geeft blijft de waarde ontbreken en maakt YieldReduction/Result daar via de afkapping min_elem(..., 1f) een volledige derving van, zodat het gewas daar niet kan landen. Daar komt ook de verziltingsopslag binnen; zie Landbouw methode.

Parameters en aannames

gxg_sign = 1: in de GxG-kaarten betekent positief een diepte onder maaiveld, dus het teken hoeft niet om; water boven maaiveld staat er negatief in.

gxg_scale = 100: de GxG-kaarten staan in meters en worden naar cm omgezet, de eenheid van de assen van de relatiedatabase. De opgelegde standen van een veenbouwsteen staan in dezelfde eenheid en met hetzelfde teken.

Zichtjaar_klimaat volgt ModelParameters/Landbouw/Scenario: de referentieband voor de kaart van 2017, de scenarioband voor 2050 en 2085.

dmgtot is schade, niet relatieve opbrengst: 0 is gezond, en de derving is de waarde gedeeld door 100.

Varianten en uitvoer

De kaartenset wordt uitgevraagd met combine_uint8 over drie assen (OpbrengstendervingsVarianten): de hydrologische leveringen NL2120_Varianten, de ankerjaren WWL_Zichtjaren en de 9 WWL_Gewassen. NL2120_Varianten is ondanks zijn naam geen variantenlijst maar de lijst met leveringen; welke variant welke leest staat in VariantK. In de NL2120-toepassing zijn dat drie leveringen, dus 3 maal 3 maal 9 is 81 kaarten. Beregening is geen as, die zit per cel in de sleutel.

De template OpbrengstenDerving_T wordt via for_each_ne voor elke combinatie geinstantieerd (container OpbrengstenDerving_Dynamisch, padstructuur <levering>/<ankerjaar>/<gewas>); de GxG-paden volgen uit de leverings- en ankerjaarnaam.

De uitvoer wordt per combinatie weggeschreven als GeoTIFF:

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

Het pad staat onder %LocalDataProjDir% en niet in de gedeelde brondata, want dit is modeluitvoer en geen brondata. 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 (Read_Opbrengstderving) en niet het rekenende item (Write_Opbrengstderving); de schakelaar is ModelParameters/OpbrengstdervingOntkoppeld. Er is geen automatische invalidatie: na een wijziging in de grondwaterstanden, de BOFEK-kaart, de beregeningskaart of de relatiedatabase moeten de kaarten opnieuw gemaakt worden, via de haak OpbrengstenDerving_Dynamisch/generate of per levering via <levering>/Generate. Die haken roepen met opzet de schrijfkant aan en niet de leeskant; WriteVariantData roept ze per variant aan.

Bekende beperkingen

De GxG en de beregeningskaart staan op 250 m, dus de grondwaterkant van de derving is blokkerig op 250 m, terwijl de bodem op 25 m varieert.

De GxG-levering dekt de Waddeneilanden niet en mist aan de randen, samen 14.600 ha landbouw van het basisjaar, en BOFEK heeft geen code op 35.000 ha landbouw. Beide worden gevuld met de dichtstbijzijnde cel die wel een waarde heeft; voor de eilanden komt die stand over zee van de Kop van Noord-Holland en de Friese kust, en dat is een tussenoplossing totdat de hydrologische levering de eilanden dekt. Een cel die ondanks de vulling geen waarde krijgt zou in de landbouwgeschiktheid een volledige derving krijgen, zodat de WWL-gewassen daar niet kunnen landen; met de vulling is dat geen enkele landbouwcel van het basisjaar.

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

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

Verantwoording

De keuzes hierboven zijn gemaakt in de volgende issues; daar staat de afweging en de meting die eraan ten grondslag ligt.

• #263: de bandindex van de relatiedatabase eengebaseerd, de opbrengstderving lineair tussen de WWL-ankerjaren, en een levering die de vernatting van de eigen variant meeneemt.
• #500: de opzet van de dynamische berekening uit de relatiedatabase, met de beregening per cel in de sleutel.
• #539: opbrengstdervingskaarten per gewas als invoer van de landbouwgeschiktheid.
• #542: beregening volgens de kaart Basisprognoses 2018 als standaardstand.
• #553: de kaarten als modeluitvoer onder %LocalDataProjDir% met de studiegebiedsuffix, en de schakelaar OpbrengstdervingOntkoppeld tussen schrijf- en leeskant.
• #780: het beregeningsregime volgt de parameter BeregeningScenario_name, en de hydrologische levering per variant in VariantK.
• #782: de klimaatband van de schaderelatie als scenariokenmerk, het basisjaar op de referentieband, de veenoplegging in de grondwaterstanden, en de open vraag over een variant op referentiestanden.
• #793: geen beregening op de cellen waar de veenbouwsteen vernat.
• #837: de GLG-as van de opzoeking ongespiegeld, de gestapelde banden via de permutatie van de opslagvolgorde bij het juiste punt gelegd, en de cellen zonder BOFEK-code of zonder grondwaterstand gevuld met de dichtstbijzijnde cel; de dekking van de eilanden door de hydrologische levering ligt bij Deltares.