diff --git a/data/fields_jena.tif b/data/fields_jena.tif index 839afe968b6e110aab2f8ea9dd44274c01ace912..1d5e88a54c96f855db61aa135e1f6c62dd486f26 100644 Binary files a/data/fields_jena.tif and b/data/fields_jena.tif differ diff --git a/data/fields_thueringer_becken.tif b/data/fields_thueringer_becken.tif new file mode 100644 index 0000000000000000000000000000000000000000..13215a12706d86af059dd631404182e54d94ad18 Binary files /dev/null and b/data/fields_thueringer_becken.tif differ diff --git a/docs/gis.md b/docs/gis.md new file mode 100644 index 0000000000000000000000000000000000000000..52ef3a804ecaa146021ae5f5a35f95d68f62c053 --- /dev/null +++ b/docs/gis.md @@ -0,0 +1,72 @@ +# GIS data in Persephone + +Persephone currently requires two separate map input files: one for land cover, +the other for field geometry. This documents describe how to obtain and process +the data needed for each of these. + +## Land cover maps + +Land cover maps for Germany at 10m resolution can be obtained from +[Mundialis](https://data.mundialis.de/geonetwork/srv/eng/catalog.search#/metadata/9246503f-6adf-460b-a31e-73a649182d07). +These are generated annually from Sentinel data and comprise the following +land cover classes: + +``` +10: forest +20: low vegetation +30: water +40: built-up +50: bare soil +60: agriculture +``` + +To create a Persephone map input file, you need to crop the national Mundialis +map to the extent that you want to simulate (suggestion: approx. 10x10km is a +reasonable size). + +To do so, download the Mundialis map and import it into QGIS. Then create a new +vector layer and create a rectangle feature to delimit the extent of your +region. Then go to `Raster -> Extraction -> Clip Raster by Extent`. Select +the Mundialis map as the input layer, set the clipping extent by choosing your +region vector layer under `Calculate from Layer` and specify the output +file name before clicking `Run`. This will generate a TIF file that you can +pass to Persephone as the `landcovermap` parameter. + +## Field ID maps + +In addition to the land cover data explained above, Persephone also needs +information about agricultural field boundaries in order to assign these to the +farming agents. Unfortunately, getting this is rather more complicated. + +In the EU, every country runs a Land Parcel Information System (LPIS) to +administer CAP payments. In Germany, this is called InVeKoS and is run by the +Länder. For example, you can view and download the InVeKoS data for Thüringen +[here](https://thueringenviewer.thueringen.de/thviewer/invekos.html). +This gives you a vector layer which can be loaded into QGIS. However, it needs +to be converted to a raster layer and cropped to your region extent before it +can be used in Persephone. + +The first thing to do is to make sure that the vector layer has a numeric (!) +field with a unique identifier for each field block (check the attribute table). +The Thüringen data has the FBI ("Feldblockident") field, but this is a string +value and therefore not usable by the rasteriser. So, we set the vector layer to +edit mode, open the field calculator, enter the information for a new field +(call it "FID" and set it to a 32-bit integer), and enter `@row_number` in the +expression field. Then save the layer and close the calculator. + +Secondly, you need to filter out all non-field/non-grassland plot types. (LPIS +also has data on forests and various landscape elements that are not relevant to +our use case.) Assuming you're working with the Thüringen InVeKoS data (other +data sets may have a different structure), right-click on the layer name in +QGIS' layer overview and click on "Filter...". Then, enter this expression in +the query builder: `"BNK" = 'AL' OR "BNK" = 'GL'` and click "OK". This +will select only field and grassland plots. + +Next, open the rasteriser (`Raster -> Conversion -> Rasterize`). Select your +FID field as the "Field to use for a burn-in value", and your land cover map (as +created above - this ensures the two layers match) as the output extent. Make sure +the "fixed value to burn" is "Not set". Then choose "Georeferenced units" as the +"Out raster size units" and set horizontal and vertical resolution to 10.0. In +the advanced parameters, set the output data type to `UInt32`. Finally, enter an +output file name and run. +