---
title: "Replace land cover and land use information in GLOBIOM"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Replace land cover and land use information in GLOBIOM}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


After creating the updated GLOBIOM input files, several adjustments need to be made in GLOBIOM to replace the country-level land cover and crop distribution information. 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: Download 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](http://www.globiom.org/) 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](https://help.github.com/articles/cloning-a-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](https://github.com/iiasa/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](https://github.com/iiasa/GLOBIOM_FABLE_Data/blob/master/README.md) 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](https://github.com/iiasa/GLOBIOM_GUI/) repository. For more information see the [GLOBIOM GUI Overview](https://iiasa.github.io/GLOBIOM_FABLE/GUI.html) 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](https://iiasa.github.io/GLOBIOM_FABLE) 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](https://help.github.com/articles/about-branches/) 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](http://www.globiom.org/) pre-compilation scripts and raw data that can be used to prepare the `Model/finaldata` input files for [GLOBIOM_FABLE](https://github.com/iiasa/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](mailto:globiom.fable@iiasa.ac.at).

#### Getting Started

To make this repository work with [GLOBIOM_FABLE](https://github.com/iiasa/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 the``Data`` 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](https://github.com/iiasa/GLOBIOM_FABLE_Data/commits/master) and [GLOBIOM_FABLE](https://github.com/iiasa/GLOBIOM_FABLE/commits/master) should be used together. For the `master` branches, compatible commits have been marked with `match_*` release tags.

As of release [match_4](https://github.com/iiasa/GLOBIOM_FABLE_Data/releases/tag/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](https://www.r-project.org/).
* The [tidyverse](https://www.tidyverse.org/) R package collection.
* The [gdxrrw](https://www.gams.com/latest/docs/T_GDXRRW.html) 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 `landcover`which 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](http://www.fao.org/3/a-i3765e.pdf).

__Table 1: USDA fertilizer farm prices__

Fertilizer |farm price [USDA-ERS](http://www.ers.usda.gov/Data/FertilizerUse/)
-----------|---------------------
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.