Replacing land cover and land use information in GLOBIOM

With respect to the distribution of crops, GLOBIOM uses two related data layers. A table with area information on the land cover per simu, distinguishing seven major land use classes,some of which are split into subclasses:

  • Cropland (CrpLnd)
  • Other agricultural land (OthAgri)
  • Grassland (Grass)
  • Forest (Forest), separated into managed, unmanaged and rotation forest
  • Wetlands (WetLnd)
  • Other natural vegetation (OthNatLnd)
  • Not relevant land cover (NotRel)

The second table provides information on the location of the 18 crops that are modeled by GLOBIOM, which together make up the Cropland class. Other agricultural land shows the location of all the other crops. The distribution of crops in the global version of the model is based on the global Spatial Production Allocation Model (SPAM) for 2000.

The function create_globiom_input() creates two gdx files. One file updates the global land cover data in GLOBIOM for the target country, while the other file replaces global land use data. Both files will be saved in the processed_data/results folder. Before create_globiom_input() can be run, you need to prepare three input files.

First, you need to collect a new country level land cover map. Any product can be used as long as it contains information on the six GLOBIOM land cover classes. The most obvious choice would be to take a national land cover map for the year 2000 or any other map that is close to the year for which the subnational statistics are available. If such map is not available it is also possible to use a global land cover product and use the country polygon to mask the relevant area. We illustrate this case for Malawi, where we used the ESACCI land cover map that was described in the section on Input data collection. We started by clipping the country area by running the script 03_spatial_data/select_esa.r (see Process spatial data).

Second, we prepared a table (mappings/esacci2globiom.csv) that maps the ESACCI land cover classes to the six GLOBIOM land cover classes. This table will be used to aggregate the land cover classes in the new map to the GLOBIOM land cover classes. You will have to prepare a similar table for your land cover map.

Finally, in case you want to add an additional crop, you have to update the table that maps the SPAMc crops to the 18 GLOBIOM crops (mappings/crop2globiom.cvs). Also see the Adding a new crop to GLOBIOM. In the example below, we will add coffee as a separate crop. The updated table will be loaded by create_globiom_input() to aggregate the SPAMc crops to the GLOBIOM crop classification.

# Select the new land cover map, in this case ESACCI
lc_file <- file.path(param$spam_path,
  glue("processed_data/maps/cropland/{param$res}/esa_raw_{param$year}_{param$iso3c}.tif"))
# Load the map
lc_map <- raster(lc_file)
plot(lc_map)

# Update crop2globiom.csv and map coff to coff (instead of the rest category) and overwrite
load_data("crop2globiom", param)
crop2globiom <- crop2globiom %>%
  mutate(globiom_crop = ifelse(crop == "coff", "Coff", globiom_crop))
write_csv(crop2globiom, file.path(param$spam_path, "mappings/crop2globiom.csv"))

# Load mapping of lc classes to globiom lc classes
lc_class2globiom <- read_csv(file.path(param$spam_path, "mappings/esacci2globiom.csv"))

# Aggregate land cover map to GLOBIOM land cover classes at simu level
# Not that the area will be expressed in 1000 ha, which is common in GLOBIOM!
create_globiom_input(lc_class2globiom, lc_map, param)

The existing land cover data in GLOBIOM can be found in the parameter LANDCOVER_INIT which defines the land cover classes at the simulation unit level and the existing crop specific land use data found in CROP_DATA parameter at the simulation unit level.

The steps required to update the GLOBIOM land cover and land use maps with new maps are the following:

  1. Download the GLOBIOM model and GLOBIOM data branches

  2. Set the regional definitions (myRegion) and the spatial resolution (0.5 degree grid, simulation unit, or LUID) and run the data compliation (0_execute_total.gms) from the data branch

  3. Start a new GAMS script that will load the appropriate files containing region and item set definitions and other national level crop production and price statistics (e.g. FAOSTAT data) from the finaldata subfolder of the model folder

  4. Load the existing data parameters for cropland use (CROP_DATA), land cover (LANDCOVER_INIT_CrpGrsForOagADJ), and load the gdx outputs from mapspam2globiom: the new cropland use crop_area and new land cover landcover maps

  5. Replace the existing cropland use in the CROP_DATA parameter with the new dataset and replace the land cover data in LANDCOVER_INIT_CrpGrsForOagADJ

  6. Identify and fill gaps in CROP_DATA for missing data for yield, water, and nutrient requirements

  7. Replace production cost data

  8. Calibrate and harmonize the cropland area and harmonize land cover area with national statistics

  9. Calibrate and rescale the crop yields using the FAO statistics

  10. Write the new parameters to the finaldata subfolder of the GLOBIOM model branch

Step 1: Downlad the GLOBIOM model and GLOBIOM data branches

For guidance in how to download the GLOBIOM model branch used for the FABLE training: https://github.com/iiasa/GLOBIOM_FABLE/blob/master/README.md

GLOBIOM_FABLE

This repository holds the GLOBIOM model code released for the FABLE training.

Getting Started

During the FABLE training, the model Model/finaldata input files were distributed via DropBox as the GLOBIOM_finaldata.zip archive. That same data has now been included in this repository.

To run the model, clone the repository to your local machine. The stages of GLOBIOM can then be executed using Model\0_executebatch.gms.

Pre compilation

The above-mentioned finaldata is mostly output from a pre-compilation stage. Pre-compilation turns raw source data into model-ready input data (finaldata). A separate GLOBIOM_FABLE_Data repository has been provided that contains the source data and pre-compilation scripts. It is required to make involved changes such as isolating a country as a separate GLOBIOM region.

To request access to the GLOBIOM_FABLE_Data repository, please email to globiom.fable@iiasa.ac.at. On acceptance of your request, you will receive an invitation from GitHub to join the GLOBIOM_FABLE_Data repository. Then you can read these instructions on how to use the data. Note that this link will be accessible only after having received access, and when signed in to GitHub.

Graphical User Interface

An optional full-featured graphical user interface (GUI) for GLOBIOM is available from the GLOBIOM_GUI repository. For more information see the GLOBIOM GUI Overview page.

Documentation

GLOBIOM is a large model comprised of many GAMS source files. To help you find your way around the model, a separate documentation website has been provided. Both the model as well as the GAMS code are documented there.

Platforms

The model has been tested to run on Windows, Linux, and MacOS. Beware that the solver can yield results that differ between platforms even when using the same GAMS version.

Terms of Use

This repository contains the GLOBIOM model as provided for the purposes of the FABLE project. This version will be regularly updated. GLOBIOM will be released under an Open Source license later. Therefore FABLE partners are requested to not redistribute the GLOBIOM model contained in this repository.

Cooperation with developing GLOBIOM is possible via this repository. Please contact us, and you will receive a repository branch with write access where you can commence with sharing your modifications and developments.

For guidance in how to download the GLOBIOM data branch:

GLOBIOM_FABLE_Data

This repository holds the GLOBIOM pre-compilation scripts and raw data that can be used to prepare the Model/finaldata input files for GLOBIOM_FABLE.

Terms of Use

This is not a public release. This repository is made available to early-adopters on an as-needed basis, and should not be shared nor redistributed: licensing terms and intellectual property rights of the various data sets are being worked out.

The content of this repository is intended for further improvement of country-specific GLOBIOM versions. If you intend to use this data for other purposes, please consult with the GLOBIOM team at IIASA.

Getting Started

To make this repository work with GLOBIOM_FABLE, clone it to a Data directory co-located with the GLOBIOM Model directory. From the command line, this can be done as follows:

cd <directory where you cloned GLOBIOM_FABLE>
git clone https://github.com/iiasa/GLOBIOM_FABLE_Data Data

This should yield a Data directory adjacent to the Model directory. Next, enter the Data directory and switch to the branch of your country team:

cd Data
git checkout <name of your branch>

Note that this will cause this Data repository to be nested inside the GLOBIOM_FABLE repository on your local machine. To make this work without mix-ups, make sure that theData directory is being ignored by the root-level .gitignore file of the GLOBIOM_FABLE repository.

To run the pre-compilation, execute Data/0_executebatch_total.gms. Pre-compilation will overwrite a subset of the files in the Model/finaldata folder.

Requirements

The pre-compilation and the model are interdependent, and as such compatible commits in GLOBIOM_FABLE_Data and GLOBIOM_FABLE should be used together. For the master branches, compatible commits have been marked with match_* release tags.

As of release match_4, pre-compilation works also on non-Window platforms and non-Latin locales. This was achieved by means of a custom developed xl2gdx.R R script that requires: * A recent version of R. * The tidyverse R package collection. * The gdxrrw R package.

Check that the Rscript binary can be invoked from the command line, and if not adjust your PATH environment variable to point to where the R binaries can be found.

Step 2: Set the regional definitions and the spatial resolution run the data compliation of the data branch

For the second step, the user should clarify which spatial resolution which will be used: the 0.5 degree grid (often called the COLROW30 resolution) or the simulation unit resolution. The simulation unit resolution (often called simu), uses soil, slope, and altitude class maps which are intersected on the 0.5 degree grid resulting in the highest spatial resolution possible in GLOBIOM.

Spatial resolution

The spatial resolution should be set and read in from the decl_Rset.gms file. Depending on the resolution of set by the user’s region/country of choice, the land use and land cover maps produced by mapspam2globiom should be aggregated to the 0.5 degree grid or left at the simulation unit level accordingly. The user should pay particular attention to make sure the correct country is selected (see code below for an example for Malawi) and that the spatial resolution is set properly. To use the Rsets.gms for a simulation unit resolution effectively some changes to this file may be necessary.

Changing the decl_Rsets.gms file for simulation unit resolution (example: Malawi)

* Define regional resolution
$setGlobal REGION REGION37
* Alternatives: GGIREGION(11) or REGION28(28) or REGION30(30) or REGIONADB(31)

$ifThen.A not setGlobal myREGION
$  ifThen.B exist "%system.FP%temp%system.dirSep%GUI_region_settings.gms"
*    Use settings from GUI
$    include "%system.FP%temp%system.dirSep%GUI_region_settings.gms"
$  else.B
*    Set country/region to single out
$    setGlobal myREGION MalawiReg
*    Use SIMU (30X30 acrmin and HRU classification) resolution (SIMU)
*    Use CR30 (30x30 arcmin ColRow) resolution (CR30)
*    USE LUID (2x2 deg) resolution (LUID) for myREGION
$    setGlobal resREGION SIMU
$  endIf.B
$endIf.A

*IMPORTANT: in the rest of the code in this file, replace the old global switch "crREGION" with new name "resRegion"


...
*Other part of the file unchanged
...

SET
SIMU_COUNTRY(COUNTRY)
CR30_COUNTRY(COUNTRY)
LUID_COUNTRY(COUNTRY)
HRUN_COUNTRY(COUNTRY)
;

*for colrow resolution
$ifthen %resREGION% == CR30
SIMU_COUNTRY(COUNTRY)
 = NO;

CR30_COUNTRY(COUNTRY)
 $ REGION_MAP('%myREGION%',COUNTRY)
 = YES;

LUID_COUNTRY(COUNTRY)
 $ SUM(REGION_MAP(REGION,COUNTRY)
         $(NOT(SAMEAS(REGION,'%myREGION%'))),1)
 = YES;
$endif

*for LUID resolution
$ifthen %resREGION% == LUID
SIMU_COUNTRY(COUNTRY)
 = NO;

CR30_COUNTRY(COUNTRY)
 = NO;

LUID_COUNTRY(COUNTRY)
 $ SUM(REGION_MAP(REGION,COUNTRY),1)
 = YES;
$endif

*for simulation unit resolution
$ifthen %resREGION% == SIMU
SIMU_COUNTRY(COUNTRY)
 $ REGION_MAP('%myREGION%',COUNTRY)
 = YES;

CR30_COUNTRY(COUNTRY)
 = NO;

LUID_COUNTRY(COUNTRY)
 $ SUM(REGION_MAP(REGION,COUNTRY)
         $(NOT(SAMEAS(REGION,'%myREGION%'))),1)
 = YES;
$endif

After the regional defintion is set, the user should then run the entire data compilation (0_execute_total.gms) in the data folder before attempting to navigate through the rest of this guide. Updating the landcover and land use maps is the last step of the data compliation before moving on to run the model. The datasets that are needed to update the cropland use and land cover maps are aggregated at the proper resolution by the data compilation file (0_execute_total.gms).

Step 3: start a GAMS script (.gms) to load the set and regional definitions and other data from the finaldata subfolder of the model folder

Start a new GAMS .gms file in the main folder of the data branch that will load the appropriate files containing region and item set definitions and other national level crop production and price statistics (e.g. FAOSTAT data) from the finaldata subfolder of the model folder

Other files should be included which contain the additional sets, parameters ,and datasets needed to properly include the new land cover and land use maps. Many of these can be found in the model folder or the finaldata folder within the model folder which will have been updated based on the data compilation performed in Step 2.

The following files should be read in by using the $include in the GAMS script:

decl_sets.gms contains all the major set definitions

decl_regionset.gms contains the regional set definitions

decl_Rsets.gms contains the switches to aggregate the model to different spatial resolutions (see Spatial Resolution section above)

sets_colrow.gms contains important spatial colrow and simulation unit information data

Load national level crop production and price statistics

This guide uses the FAOSTAT crop production data from proddata_c.gms and the FAOSTAT crop price data from data_supply.gms as its main harmonization dataset however any good national level crop production and price statistic dataset with a full information on crop area, crop production, and price can be used. Load into the .gms file the two FAOSTAT files found in the finaldata subfolder in the model folder. PRODDATA_C is the name of the parameter which contains the national level production statistics from FAO which will be used for rescaling the crop yields and SData is the parameter which contains the national level crop price data.

Example code for including GAMS (.gms) files

*load the set definitions 
$include decl_sets.gms 

Step 4: Loading and replacing the cropland use maps in the appropriate parameters within GLOBIOM

Load the existing cropland use maps

In the GAMS .gms file, load decl_paramgdx.gms and data_crops.gdx which are found in the finaldata subfolder in the model folder.

decl_paramgdx.gms properly defines the CROP_DATA parameter. The CROP_DATA parameter and its associated data is stored in the .gdx file.

The CROP_DATA parameter is the parameter which stores the crop land use are and production management information. Updating the CROP_DATA parameter with the new crop land use maps and other crop production data is the first objective of this guide.

Load the existing land cover maps

In the GAMS .gms file, load the land cover parameter found in the data_EPICLandCov_CrpGrsForOagADJ.gms which is found in the finaldata subfolder in the model folder and is compiled by the data compilation using the existing datasets.

The parameter which contains the existing landcover data by land unit is called LANDCOVER_INIT_CrpGrsForOagADJ. Updating the LANDCOVER_INIT_CrpGrsForOagADJ parameter with the new land cover maps is the second objective of this guide.

Load the new cropland use map and land cover maps from mapspam2globiom

In the GAMS .gms file, load the new landcover dataset found in globiom_landcover_YEAR_ISO3.gdx(YEAR and ISO3 are replaced by the user’s ISO3 country code and the year of data)

This dataset is produced as an output from the mapspam2globiom R package. The only parameter in this .gdx file is a parameter called landcoverwhich provides the new area information on the land cover area per simulation unit is distinguished by six different land use categories:

• Total crop and other agricultural land (CrpLnd, OthAgri)

• Grassland (Grass)

• Forest (Forest)

• Wetlands (WetLnd)

• Other natural vegetation (OthNatLnd)

• Not relevant land cover (NotRel)

Together with the grassland class (Grass), the total crop and other agricultural land classes comprise the total agricultural land cover areas by simulation unit within GLOBIOM. The land cover class total crop and other agricultural area is divided into two parts: (1) Cropland (CrpLnd), which presents the location of the 20 crops that are modelled by GLOBIOM and (2) Other agricultural land (OthAgri), which presents the location of all the other crops not explicitly modeled by GLOBIOM.

The total crop and other agricultural land class is linked with the data layer with a more detailed information on the cropland use by crop and management system per simu called globiom_crop_area_YEAR_ISO3.gdx (year and ISO3 are replaced by the user’s ISO3 country code and the year of data) produced as an output from the mapspam2globiom R package. In this .gdx file there is a parameter called crop_area which is the new cropland use area dataset.

Load globiom_crop_area_YEAR_ISO3.gdx into the GAMS file.

Example code for including gdx files

*define the parameter 
parameter
crop_area(ANYREGION,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,ALLTECH) new land use map from mapspam2globiom
;

*load new cropland use map from mapspam2globiom package
$GDXIN globiom_crop_area_YEAR_ISO3.gdx
$LOADDC crop_area
$GDXIN

Step 5: Replace existing crop land use in CROP_DATA and replace the existing land cover data in LANDCOVER_INIT_CrpGrsForOagADJ

Add code to the GAM .gms file that will save the existing CROP_DATA and LANDCOVER_INIT_CrpGrsForOagADJ information for the country into a separate parameter to be used later for harmonization wit the national level statistics.

Parameter
CROP_DATA_Orig;

CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,ALLTECH,ALLITEM)  =
CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,ALLTECH,ALLITEM);
*Note here: SIM_COUNTRY is used for a simulation unit resolution 

Parameter
LANDCOVER_Orig;

LANDCOVER_Orig(SIMU_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC)= 
LANDCOVER_INIT_CrpGrsForOagADJ(SIMU_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC); 

Add to the code of the GAMS .gms file, replace the existing CROP_DATA cropland use area for the country with the new cropland use dataset (crop_area) (example for the simulation unit resolution below).

CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,INPUTSYS,'BaseArea') 
= crop_area(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,INPUTSYS);

Add to the code of the GAMS .gms file, replace the existing CROP_DATA cropland use area for the country with the new cropland use dataset (crop_area). An example for COLROW resolution is below which replaces the existing CROP_DATA, with new dataset while also aggregating to the COLROW/CR30 resolution.

parameter
crop_area2(ANYREGION,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,ALLTECH) CR resolution of the new crop areas
; 

crop_area2(CR30_COUNTRY,COLROW30,"Alti_Any","Slp_Any","Soil_Any",AEZCLASS,SPECIES,INPUTSYS)= sum((ALTICLASS,SLPCLASS,SOILCLASS),croparea_spam(CR30_COUNTRY,COLROW30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,INPUTSYS));
*Note here: CR30_COUNTRY is used for the 0.5 degree COLROW resolution 

CROP_DATA(CR30_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,INPUTSYS,'BaseArea') = crop_area(CR30_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,INPUTSYS);

For simulation unit resolution, add code that will replace the existing LANDCOVER_INIT_CrpGrsForOagADJ land cover area for the country with the new land cover use dataset (land_cover).

LANDCOVER_INIT_CrpGrsForOagADJ(SIMU_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC) 
= landcover(SIMU_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC);

For COLROW resolution, add code to the .gms file that will replace the existing CROP_DATA, with new dataset while also aggregating to the COLROW/CR30 resolution

parameter
landcover2(CR30_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC) CR resolution of the new landcover
; 

landcover2(CR30_COUNTRY,COLROW30,"Alti_Any","Slp_Any","Soil_Any",AEZCLASS,SPECIES,INPUTSYS) = sum((ALTICLASS,SLPCLASS,SOILCLASS), croparea_spam(CR30_COUNTRY,COLROW30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,INPUTSYS));

*Note here: CR30_COUNTRY is used for the 0.5 degree COLROW resolution 

LANDCOVER_INIT_CrpGrsForOagADJ(CR30_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC)
= landcover2(CR30_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC);

Troubleshooting: reorder or drop dimenstions

If the mapspam2globiom output has too many dimensions to fit in CROP_DATA or LANDCOVER_INIT_CrpGrsForOagADJ or the dimensions are in the wrong order you can reorder a parameter like this:

parameter
crop_area_reorder;

*put new croparea spam data in correct order
crop_area_reorder(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,INPUTSYS)
         $LANDCOVER_INIT_CrpGrsForOagADJ(SIMU_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,'SimUarea')=
crop_area(SPECIES,INPUTSYS,SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS);

parameter
crop_area_dropdim;


*drop the simulation unit number from the new crop_area spam data
crop_area_dropdim(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,SPECIES,INPUTSYS)
         $LANDCOVER_INIT_CrpGrsForOagADJ(SIMU_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,'SimUarea')=
sum(SimUID, crop_area(SimUID,SPECIES,INPUTSYS,SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS));

Step 6: Udpate the crop yields for new crop areas and identify and fill missing data with averages values

To use the cropland use area map, cropland management and crop production information is needed in GLOBIOM to model the production at grid cell or unit level. A dataset of crop yields and input requirements should provide information that can be linked to the crop/management systems in each simu location. The dataset of yields and input requirements can be produced using outputs of regional biophysical crop models or locally reported yield and input data. The yield and input requirements dataset can also be generated using the existing crop model outputs available for GLOBIOM. GLOBIOM is likely to have existing biophysical crop model yields and input requirements for many crop/management systems at a grid cell or simulation unit level even if the previous cropland use map does not show crop area in that gridcell or simulation unit location in the base year. However if the newly generated cropland use map shows new crops/management systems in simu locations with no existing crop model yields and input requirements, it is necessary to fill this information gap. To do this we replace the missing information with the average national average values.

Add code to replace the CROP_DATA crop yields for the new crop areas:

CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,PRODUCT)=
CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,PRODUCT);

Calculate the crop yield average using FAO Statistics

parameter
average_FAO_yield(anyregion,crops)
;
*calculate the average yield based on the FAO statistics
average_FAO_yield(simu_country,crops)$PRODATA_C(simu_country,CROPS,'AreaHarv') =
 (PRODATA_C(simu_country,CROPS,'Yield')*PRODATA_C(simu_country,CROPS,'AreaHarv'))/
 PRODATA_C(simu_country,CROPS,'AreaHarv');

Calculate the crop yield average using the average of each management system (based on EPIC Data)

parameter
average_EPIC_yield(anyregion,inputsys,crop)

;
*calculate the average yield for each system based on the EPIC data and existing land use map
average_EPIC_yield(simu_country,inputsys,crop) $ sum((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS), CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,'BaseArea')) = sum((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROPS), CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,CROPS)* (CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,'BaseArea')))/ (sum((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS), CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,'BaseArea')));

Calculate the average crop input requirement using the average of each management system (based on EPIC Data)

parameter
average_FT(anyregion,inputsys,crop,*)
;
*Note: nitrogen = 'FTN'; phosporus = 'FTP' and irrigation water requirement = 'WATER' 

*calculate the average nitrogen by crop and system
average_FT(simu_country,inputsys,crop,'FTN')$ sum((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS), CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,'BaseArea')) = sum((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS),  CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,'FTN')* (CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,'BaseArea')))/ (sum((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS), CROP_DATA_Orig(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,'BaseArea')));

Identify the grid cells with missing data and replace with average crop yields

Use the following code to identify the grid cells or units with missing crop yield data and then replace it with the country average crop yield or other data sources.

Parameter
missing_crop_data;

*identify missing crop yield data (first pass) 
missing_crop_data(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,"Problem_1")
 $(CROP_DATA(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,"BaseArea") and
 (sum(crops,CROP_DATA_Orig(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,crops))=0))=
 1;
 
*only replace the missing crop yields with management system specific crop yield averages when there is missing data after first pass 
CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,PRODUCT)$ 
 $(CROPPRODMAP(CROP,PRODUCT) and missing_crop_data(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,"Problem_1"))=
  average_EPIC_yield(SIMU_COUNTRY,inputsys,crop);
 
 
*identify missing crop yield data (second pass)
missing_crop_data(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,"Problem_2")
 $(CROP_DATA(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,"BaseArea") and
 (sum(crops,CROP_DATA_Orig(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,crops))=0))=
 1; 
 
*after second pass replace the missing crop yields with national averages
CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,PRODUCT)$ 
 $(CROPPRODMAP(CROP,PRODUCT) and missing_crop_data(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,"Problem_2"))=
  average_FAO_yield(SIMU_COUNTRY, PRODUCT);
  
*identify missing crop yield data (third pass)
missing_crop_data(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,"Problem_3")
 $(CROP_DATA(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,"BaseArea") and
 (sum(crops,CROP_DATA_Orig(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,inputsys,crops))=0))=
 1; 
 
 
 *if there are still missing crop yields after the third pass find an additional source or use a world average for crop yields. 

Step 7: Update production costs by crop and management system

Rather than rely on the exisiting dataset for crop production costs from CROP_DATA, it is eaiser to recalculate all the costs after updating the crop areas. The production costs are calculated with the assumptions that for low and subsistence crop production, the total costs equal the total revenues per ha. For high input, rainfed we also add the additional nutrient expenses based on USDA fertilizer farm prices, based on the Table 1. For irrigation systems, we add the additional nutrient requirement and also operations and maintenance costs for irrigation systems based on an FAO estimate based on smallholder irrigation Smith et al., 2014.

Table 1: USDA fertilizer farm prices

Fertilizer farm price USDA-ERS
Nitrogen 0.589 USD per kg
Phosphorus 0.562 USD per kg

Table 2

Irrigation Systems Operations and maintenance costs
Basin 370 USD per ha
Furrow 370 USD per ha
Sprinkler systems 1200 USD per ha
Drip systems 1760 USD per ha

Example code for updating the crop/management system production costs

The following code should be added to the GAMS .gms file to replace the crop production costs.

PARAMETER
PRICE_FT(FT) farm prices of fertilisers in USD per kg
* SOURCE: http://www.ers.usda.gov/Data/FertilizerUse/ ,     Table 7
* calculated as average over 2001-2005 and recalculated on pure nutrients
/
 FTN    0.589
 FTP    0.562
/

* subsistence management system costs
CROP_DATA(simu_country,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'SS','COST')
 $(CROP_DATA(simu_country,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'SS',"BaseArea") and
 sum(crops$cropprodmap(crop,crops), sum(region$region_map(region,simu_country),PRODATA_C(REGION,CROPS,'Yield'))))
 = sum(crops$cropprodmap(crop,crops), sum(region$region_map(region,simu_country),
         SData(REGION,CROPS,"Price") * PRODATA_C(REGION,CROPS,'Yield')));

CROP_DATA(simu_country,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'SS','COST')
 $(CROP_DATA(simu_country,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'SS',"BaseArea") and
 not(sum(crops$cropprodmap(crop,crops), sum(region$region_map(region,simu_country),PRODATA_C(REGION,CROPS,'Yield')))))
 = sum(crops$cropprodmap(crop,crops), sum(region$region_map(region,simu_country),
         SData(REGION,CROPS,"Price") * average_yield(crops)));

* low input rainfed management system costs
CROP_DATA(simu_country,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'LI','COST')
 $(CROP_DATA(simu_country,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'LI',"BaseArea") and
 sum(crops$cropprodmap(crop,crops), sum(region$region_map(region,simu_country),PRODATA_C(REGION,CROPS,'Yield'))))
 = sum(crops$cropprodmap(crop,crops), sum(region$region_map(region,simu_country),
         SData(REGION,CROPS,"Price") * PRODATA_C(REGION,CROPS,'Yield')));

CROP_DATA(simu_country,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'LI','COST')
 $(CROP_DATA(simu_country,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'LI',"BaseArea") and
 not(sum(crops$cropprodmap(crop,crops), sum(region$region_map(region,simu_country),PRODATA_C(REGION,CROPS,'Yield')))))
 = sum(crops$cropprodmap(crop,crops), sum(region$region_map(region,simu_country),
         SData(REGION,CROPS,"Price") * average_yield(crops)));
* high input rainfed management system costs
CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'HI','COST')
 $ CROP_DATA(simu_country,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'HI',"BaseArea")
 = sum(crops$cropprodmap(crop,crops), sum(region$region_map(region,simu_country),
            SData(REGION,CROPS,"Price") * PRODATA_C(REGION,CROPS,'Yield')))
 + SUM(FT,(CROP_DATA(SIMU_COUNTRY,COLROW30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'HI',FT)
          -CROP_DATA(SIMU_COUNTRY,COLROW30,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,'LI',FT))
  * PRICE_FT(FT) * 2.5);

Step 8: Harmonize the cropland area with national level statistics and the land cover with the national level statistics

The purpose of this step is to harmonize the new base year crop area with the reported FAO statistics for crop area. The user can use any national crop area statistics but should take care to compare only physical area as the cropland use maps produced by mapspam2globiom are in physical areas not harvested area. The following code should be added to the .gms file and used to calculate the differences in the crop areas with the national statistics.


SET
SOURCE
/
OBS_FAO  Observed data from FAOSTAT
OBS_STAT Observed data from updated SPAM dataset
OBS_SIMU  Observed data from updated GLOBIOM parameter
CALC
'PCT%'   percent difference
OBS_STATINIT Observed data from the original source 
/

PARAMETER
ACR_CHECK(ANYREGION,SPECIES,SOURCE) area check parameter in 1000 ha;

*From FAO
ACR_CHECK(SIMU_COUNTRY,CROP,'OBS_FAO')
 = SUM(PRODUCT $CROPPRODMAP(CROP,PRODUCT),
                  PRODATA_C(SIMU_COUNTRY,PRODUCT,'AreaHarv'))/1000;
                  
               
*From the new crop area
ACR_CHECK(SIMU_COUNTRY,CROP,'OBS_STAT')
 = sum((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,INPUTSYS),
    crop_area(SIMU_COUNTRY,colrow30,ALTICLASS,SLPCLASS,SOILCLASS,INPUTSYS,CROP)  ;


*from new simulation unit data in CROP_DATA
ACR_CHECK(SIMU_COUNTRY,CROP,'OBS_SIMU')
 = SUM((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,InputSys)
        $ (CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,InputSys,'BaseArea')),
     CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,InputSys,'BaseArea')) ;

ACR_CHECK(SIMU_COUNTRY,CROP,'PCT%')
 $ ACR_CHECK(SIMU_COUNTRY,CROP,'OBS_STAT')
 = (ACR_CHECK(SIMU_COUNTRY,CROP,'OBS_SIMU')/ACR_CHECK(SIMU_COUNTRY,CROP,'OBS_STAT')
    -1) *100;

ACR_CHECK(SIMU_COUNTRY,CROP,'PCT%')
 $ (abs(ACR_CHECK(SIMU_COUNTRY,CROP,'PCT%')) lt 0.001)
 = 0;

The purpose of this step is to check the differences between the total cropland cover, the new cropland cover, and the aggregated cropland use. The mapspam2globiom package uses the new cropland cover dataset as an input for the cropland use maps so this step is simply to check that all crops are included properly.

PARAMETER
LndCov_CHECK(ANYREGION,SOURCE) land cover check parameter in 1000 ha;

*From the FAO statistics 
LndCov_CHECK(SIMU_COUNTRY,'OBS_FAO')
= SUM(PRODUCT, PRODATA_C(SIMU_COUNTRY,PRODUCT,'AreaHarv'))/1000;


*From the original landcover map
ACR_CHECK(SIMU_COUNTRY,CROP,'OBS_STATINIT')
 = SUM((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS),
           LANDCOVER_Orig(SIMU_COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC) 

*From the new crop area
LndCov_CHECK(SIMU_COUNTRY,'OBS_STAT')
 = sum((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS),
landcover(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,'CrpLnd'));

*from new simulation unit data in CROP_DATA
LndCov_CHECK(SIMU_COUNTRY,'OBS_SIMU')
 = SUM((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS)
   LANDCOVER_INIT_CrpGrsForOagADJ(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,'CrpLnd')) ;


LndCov_CHECK(SIMU_COUNTRY,'PCT%','1')
 $ LndCov_CHECK(SIMU_COUNTRY,'OBS_STAT','1')
 = (LndCov_CHECK(SIMU_COUNTRY,PRODUCT,'OBS_SIMU','1')/LndCov_CHECK(SIMU_COUNTRY,'OBS_STAT','1')
    -1) *100 ;

LndCov_CHECK(SIMU_COUNTRY,'PCT%','1')
 $ (abs(LndCov_CHECK(SIMU_COUNTRY,'PCT%','1')) lt 0.001)
 = 0;

Step 9: Calibrate and rescale the crop yields using the national level statistics

The purpose of this step is to harmonize the crop production in the base year that is calculated using the new base year crop areas with the production quantitites from the national level crop produciton statistics in the base year. The following code will identify and rescale the crop yields in the CROP_DATA parameter so that the calculated crop production levels in the base year (crop yield * new crop area) match the reported statistics.

PARAMETER
PRODCHECK(ANYREGION,ALLPRODUCT,SOURCE,*)
YIELDCHECK(ANYREGION,ALLPRODUCT,SOURCE,*)
;

PRODCHECK(SIMU_COUNTRY,PRODUCT,'OBS_FAO','1')
  = PRODATA_C(SIMU_COUNTRY,PRODUCT,'ProdQ') * 0.001;

PRODCHECK(SIMU_COUNTRY,PRODUCT,'OBS_SIMU','1')
 = SUM((ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,InputSys,CROPPRODMAP(CROP,PRODUCT))
       $CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,InputSys,PRODUCT)   ,
     CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,InputSys,PRODUCT)
         *CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,InputSys,'BaseArea'));


PRODCHECK(SIMU_COUNTRY,PRODUCT,'PCT%','1')
 $ PRODCHECK(SIMU_COUNTRY,PRODUCT,'OBS_FAO','1')
 = (PRODCHECK(SIMU_COUNTRY,PRODUCT,'OBS_SIMU','1')/PRODCHECK(SIMU_COUNTRY,PRODUCT,'OBS_FAO','1')
    -1) *100 ;

PRODCHECK(SIMU_COUNTRY,PRODUCT,'PCT%','1')
 $ (abs(PRODCHECK(SIMU_COUNTRY,PRODUCT,'PCT%','1')) lt 0.001)
 = 0;

LOOP(CROPPRODMAP(CROP,PRODUCT),
YIELDCHECK(SIMU_COUNTRY,PRODUCT,SOURCE,'1')
 $ ACR_CHECK(SIMU_COUNTRY,CROP,SOURCE)
 = PRODCHECK(SIMU_COUNTRY,PRODUCT,SOURCE,'1')
 / ACR_CHECK(SIMU_COUNTRY,CROP,SOURCE) ;
);

YIELDCHECK(SIMU_COUNTRY,CROPS,'PCT%','1')
 $ YIELDCHECK(SIMU_COUNTRY,CROPS,'OBS_FAO','1')
 = (YIELDCHECK(SIMU_COUNTRY,CROPS,'OBS_SIMU','1')/YIELDCHECK(SIMU_COUNTRY,CROPS,'OBS_FAO','1')
    -1) *100;

* yield adjustment to match production in the FAO statistics
CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,InputSys,PRODUCT)
 $(CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,InputSys,PRODUCT) AND
   PRODCHECK(SIMU_COUNTRY,PRODUCT,'OBS_FAO','1')           AND
   PRODCHECK(SIMU_COUNTRY,PRODUCT,'OBS_SIMU','1') )
 = CROP_DATA(SIMU_COUNTRY,ALLCOLROW,ALTICLASS,SLPCLASS,SOILCLASS,AEZCLASS,CROP,InputSys,PRODUCT) *
   PRODCHECK(SIMU_COUNTRY,PRODUCT,'OBS_FAO','1')/ PRODCHECK(SIMU_COUNTRY,PRODUCT,'OBS_SIMU','1');

Step 10: Send the updated parameters to the finaldata subfolder in the model branch

The final step in the process to add new base year cropland use and land cover maps to use in GLOBIOM is to send the parameters to the finaldata subfolder of the model branch.

This process is relatively straightforward depending on how the data and model branches are linked on the user’s machine.

The CROP_DATA parameter is held in a GDX file. Depending on how the folders are nested on the user’s machine, the code to send the parameter to the appropriate GDX in the finaldat folder can look something like this:

$setLocal X %system.dirSep%     

execute_unload '..%X%Model%X%finaldata%X%data_crops.gdx'          , CROP_DATA

To write the landcover parameter as a .gms file with the appropriate name: data_EPICLandCov_CrpGrsForOagADJ.gms the following code should be used;

FILE Resource_DataLandDet_CrpGrsForOagADJ /'..%X%Model%X%finaldata%X%data_EPICLandCov_CrpGrsForOagADJ.gms'/;
PUT  Resource_DataLandDet_CrpGrsForOagADJ;

Resource_DataLandDet_CrpGrsForOagADJ.pw = 1000;
Resource_DataLandDet_CrpGrsForOagADJ.lw = 15;
Resource_DataLandDet_CrpGrsForOagADJ.lj = 1;
Resource_DataLandDet_CrpGrsForOagADJ.nw = 15;
Resource_DataLandDet_CrpGrsForOagADJ.nd = 2;


PUT "TABLE LANDCOVER_INIT_CrpGrsForOagADJ(ALLCOUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC) Initial land cover adjusted for consistency with FAO arable land grass req and forest (1000 ha)" /;

PUT @80;
LOOP(LC_TYPES_EPIC,
 PUT LC_TYPES_EPIC.TL; );
PUT /;

Resource_DataLandDet_CrpGrsForOagADJ.lj = 2;
Resource_DataLandDet_CrpGrsForOagADJ.lw = 0;

LOOP((COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass)
 $ SUM(LC_TYPES_EPIC,LANDCOVER_INIT_CrpGrsForOagADJ(COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC)),
 PUT COUNTRY.TL,'.',ALLCOLROW.TL,'.',AltiClass.TL,'.',SlpClass.TL,'.',SoilClass.TL,'.',AezClass.TL;
 PUT @80
 LOOP(LC_TYPES_EPIC,
  IF(LANDCOVER_INIT_CrpGrsForOagADJ(COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC) ne 0,
   PUT (LANDCOVER_INIT_CrpGrsForOagADJ(COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC)););
  IF(LANDCOVER_INIT_CrpGrsForOagADJ(COUNTRY,ALLCOLROW,AltiClass,SlpClass,SoilClass,AezClass,LC_TYPES_EPIC) eq 0,
   PUT "               ";);
     );
 PUT /;
     );
PUT ";";

Once these ten steps are taken the user can navigate toward the model folder and should run the GLOBIOM model with the new land cover and cropland use maps.