In AndySouth/climcropr: Simple Crop Modelling at a Global Scale

knitr::opts_chunk$set(echo = TRUE, fig.width=7)

library(tidyverse)
library(stringr)
# library(forcats)
# library(sp) #for maps
# library(tmap) #for maps
library(knitr) #for tables with kable
library(climcropr)
library(raster)
library(rnaturalearth) #for maps
library(tmap) # for maps
library(dismo) #for ecocrop

# extent object for use in plots
ext <- extent(-180,180,-60,80)

Draft paper. Journal options :

1. Env modelling and software : they encourage short communications of < 3000 words. 4.4.
2. Agricultural and Forest Meteorology : Ramirez & other relevant recent papers. 3.9.
3. Agricultural systems : 3 most recent Agmip papers in. 2.5.

Abstract

We describe a simple, generic, mechanistic approach for investigating the broad responses to climate change of the global distributions of hundreds of crops.

We use temperature and rainfall requirements stored in the FAO ecocrop database to assess crop suitability by comparison with predicted monthly climate data. Our method is based on, and simpler than, existing Ecocrop modelling approaches.

Our method predicts if an area is suitable or not (0/1), rather than predicting a suitability value between 0 and 1

We apply this simpler method to X crops at a global scale and test results against MIRCA harvested area data for the year 2000.

This new simpler method produces better true positive rates than the original ecocrop model when tested against the MIRCA data.

Existing ecocrop models have not been applied to a large number of crops at large geographic extents even though the potential to do this is one of the advantages of its simplicity.

This new approach is both simpler and more transparent. We provide a documented R package implementing the method, including methods to test against field data.

~ end of abstract

[@Ramirez-Villegas2013] sorghum paper [@Anderson2015] differences in global crop models and maps [@Hijmans2012] diva-gis manual [@Hijmans2017] dismo R package [@Estes2013] comparing mechanistic and empirical crop models [@Holzkamper2015] climate change impacts on maize based on 3 different modelling approaches [@Iizumi2015] how do weather influence croppping area and intensity [@Leff2004] geographic distribution of major crops [@Mueller2017] global gridded crop model evaluation [@Porter2014] IPCC AR5 agriculture chapter [@Portmann2010] MIRCA crop maps

Introduction

Changes in temperature and rainfall predicted in high-end climate scenarios are outside of the tolerances of many existing crops in many of the areas where they are currently grown ?[@Porter2014]. Potential adaptation strategies include growing crops in new areas that become suitable under a changed climate and stopping growing certain crops in areas that become unsuitable. Global assessments of likely impacts on crops and the potential for these adaptation strategies would be aided by an ability to predict suitability of areas for a wide range of crops at the global scale.

There are a wide range of existing approaches for predicting the response of crop distributions to climate change [?refs]. These differ in whether they are statistical or mechanistic and in the level of detail they include. There is a trade off between this level of detail and the generality of application in terms of numbers of crops and geographic area. Detailed models, e.g. representing the growth of crops daily, have data requirements that restrict their application to very few crops and to restricted geographic areas.

A relatively simple mechanistic approach based upon the temperature and precipitation requirements of hundreds of crops held in the FAO EcoCrop database was developed more than 15 years ago. This modeling approach, also termed EcoCrop, was implemented in both GIS and a package in the statistical language R. A more recent paper describes the application of a modified version of this approach to a single crop, Sorghum[@Ramirez-Villegas2013].

The FAO ecocrop database holds data for more than 2500 crops including requirements for temperature, rainfall and soils based upon expert knowledge. The ecocrop modelling approach estimates suitability for crops based on temperature and precipitation using these requirements[@Ramirez-Villegas2013]. The ecocrop modelling approach has been implemented in a number of slightly different ways and can be run from an R package[@Hijmans2017] and Diva-GIS[@Hijmans2012]. Ecocrop the model is designed to be simpler than other crop models allowing it to be applied to a wider range of crops and situations.

In summary the ecocrop the modelling approach : 1. Calculates a mean of the minimum and maximum crop growth cycle 2. Rounds this mean (up or down) to the nearest number of months to calculate a 'duration' 3. Finds whether there are 'duration' consecutive months in between the ecocrop temperature limits 4. Finds whether there is a concurrent 'duration' period where the total precipitation is between the ecocrop rainfall requirements of the crop.

This is a summary of the ecocrop model as implemented within the R package dismo and within Diva-GIS.

The original ecocrop model uses 'optimal' temperature and precipitation limits specified in FAO ecocrop to develop an index of suitability per month, the best month and the value of the index in the best month.

Ecocrop the model as implemented in this way has not been well tested against empirical data.

The existing ecocrop modelling approach involves a number of steps that are not entirely straightforward and, we would argue, are not well justified.

Firstly a mean of the minimum and maximum crop growth cycles is used as one of the most influential parameters. For example for maize the minimum and maximum crop growth cycles are 65 days and 365 days which round to 3 months and 12 months. The ecocrop model thus assesses temperatures for a 7 month growth cycle. The figure below shows the minimum and maximum crop growth cycles stored within ecocrop for some important crops.

Secondly the crop parameters as specified in the FAO ecocrop database are necessarily very broad. Whilst they are based on expert judgement they cannot hope to capture the majority of the processes that influence the suitability of an area for a crop. For this reason a more realistic goal than trying to predict a suitability score is to try to predict which areas are definitely not suitable for a crop. This also has the advantage that predicted absence is easier than a suitability score to test against global crop data.

We aim to move towards a simpler more straightforward method.

Methods

The proposed method is :

1. Use the minimum crop growth cycle to assess whether there is a minimum period within the year within the crop temperature limits that could allow the crop to grow.
2. Use the maximum crop growth cycle to assess whether there is a period that could allow the crop to get sufficient water. Trials have shown that using either the maximum crop growth cycle to assess temperature or the minimum crop growth cycle to assess rainfall leads to the exclusion of unrealistic areas (e.g. excluding maize from most of Europe).

The base ecocrop predictions have some aspects that also look wrong for example they too predict that maize will not grow in Europe, the ecocrop method of using the mean growing cycle to find time periods without low temperatures excludes moderately northern latitudes.

TODO ~ I need to justify advantages over a statistical approach (statistical works on where crops are grown not on where they could be grown)

Cut from above : It aims for a conservative approach, to be able to exclude areas that are unsuitable with some confidence, while being less concerned that some unsuitable areas are classed as suitable. [todo I'm not sure about that sentence]

Table 1 : The new simpler method has a better true positive rate than ecocrop for rainfed crops. Comparing global predictions for the year 2000 to Mirca.

# new simpler method for predicting suitability 

#cropnames <- c('maize','wheat','rice','potato','soya bean','sugarcane')
#beware underscore in soya_bean
cropnames <- c('maize','rice','wheat','soya_bean')
#using cr_suit_simpler


#dataframe to store results to output a table
dfres <- data.frame( crop=cropnames,
                     simple_true_positive=NA,
                     ecocrop_true_positive=NA)
# dfres <- data.frame( crop=cropnames,
#                      simple_true_positive=NA,
#                      ecocrop_true_positive=NA,
#                      simple_true_negative=NA,
#                      ecocrop_true_negative=NA)


cropnum <- 0
for( crop in cropnames)
{
  cropnum <- cropnum +1

  # this uses st, a raster stack stored in the package containing crop predictions & mirca data
  obs <- subset(st,paste0(crop,'_mirca'))

  #compare suitsimp to Mirca
  pred <- subset(st,paste0(crop,'_suitsimp'))

  tpos_s <- true_pos(obs,pred)
  dfres$simple_true_positive[cropnum] <- signif(tpos_s,2)

  #dfres$simple_true_negative[cropnum] <- signif(true_neg(obs,pred),2)  

  #compare ecocrop to Mirca
  pred <- subset(st,paste0(crop,'_ecocrop'))

  tpos_e <- true_pos(obs,pred)
  dfres$ecocrop_true_positive[cropnum] <- signif(tpos_e,2) 

  #dfres$ecocrop_true_negative[cropnum] <- signif(true_neg(obs,pred),2) 

  #cat(paste(crop, 'simp:', tpos_s, 'ecocrop:', tpos_e, '\n\n' ))

}

kable(dfres)

# testing tmax method against simpsuit and ecocrop
# based on files held in memory and saved as Rds
# so eval=FALSE for now
st_tmax <- readRDS('C:\\Dropbox\\ueaHelix2017\\ecocrop\\andy_code\\st_tmax.Rds')

#cropnames <- c('maize','wheat','rice','potato','soya bean','sugarcane')
#beware underscore in soya_bean
cropnames <- c('maize','rice','wheat','soya_bean','sorghum')

#dataframe to store results to output a table
dfres <- data.frame( crop=cropnames,
                     simple_true_positive=NA,
                     tmax_true_positive=NA,                     
                     ecocrop_true_positive=NA)

cropnum <- 0
for( crop in cropnames)
{
  cropnum <- cropnum +1

  # this uses st, a raster stack stored in the package containing crop predictions & mirca data
  obs <- subset(st,paste0(crop,'_mirca'))

  #tpos suitsimp to Mirca
  pred <- subset(st,paste0(crop,'_suitsimp'))
  tpos_s <- true_pos(obs,pred)
  dfres$simple_true_positive[cropnum] <- signif(tpos_s,2)

  #tpos tmax to Mirca
  #currently stored in different place
  pred <- subset(st_tmax,paste0(crop,'_tmax'))
  tpos_s <- true_pos(obs,pred)
  dfres$tmax_true_positive[cropnum] <- signif(tpos_s,2)

  #tpos ecocrop to Mirca
  pred <- subset(st,paste0(crop,'_ecocrop'))
  tpos_e <- true_pos(obs,pred)
  dfres$ecocrop_true_positive[cropnum] <- signif(tpos_e,2) 




  #dfres$ecocrop_true_negative[cropnum] <- signif(true_neg(obs,pred),2) 

  #cat(paste(crop, 'simp:', tpos_s, 'ecocrop:', tpos_e, '\n\n' ))

}

kable(dfres)

Table 2 : The new simpler method has a slightly poorer true negative rate (specificity) than ecocrop for rainfed crops. However the true negative rate is not such a good test because we expect that other factors (e.g. such as history and economy) contribute to crops not being grown in areas that have a suitable climate.

# new simpler method for predicting suitability 

#cropnames <- c('maize','wheat','rice','potato','soya bean','sugarcane')
#beware underscore in soya_bean
cropnames <- c('maize','rice','wheat','soya_bean')
#using cr_suit_simpler


#dataframe to store results to output a table
dfres <- data.frame( crop=cropnames,
                     simple_true_negative=NA,
                     ecocrop_true_negative=NA)


cropnum <- 0
for( crop in cropnames)
{
  cropnum <- cropnum +1

  # this uses st, a raster stack stored in the package containing crop predictions & mirca data
  obs <- subset(st,paste0(crop,'_mirca'))

  #tneg suitsimp to Mirca
  pred <- subset(st,paste0(crop,'_suitsimp'))
  dfres$simple_true_negative[cropnum] <- signif(true_neg(obs,pred),2)  

  #tneg ecocrop to Mirca
  pred <- subset(st,paste0(crop,'_ecocrop'))
  dfres$ecocrop_true_negative[cropnum] <- signif(true_neg(obs,pred),2) 

}

kable(dfres)

Maps showing how the new simpler predictions compare both to Mirca and the existing ecocrop model.

Both Mirca and ecocrop have a gradation of values, these are converted to 0/1 in the testing process.

In the first set of maps I convert everything to 0/1 to aid comparisons.

In the 2nd set of maps I keep the range of values in the Mirca and ecocrop predictions. In these later maps grey in the Mirca maps represents harvested area greater than 0 (but low). Grey in the suitability maps represents 0 suitability. I think that the fact that the green pattern in the suitability maps follows the grey pattern in the mirca maps is a good sign. I just want us to discuss this and then I can improve the colour schemes.

\pagebreak

---

#fig.height reduced to try to get 3 plots on a page

  #beware underscore in soya_bean
  cropnames <- c('maize','rice','wheat','soya_bean')

  #this doesn't currently help with plotting maps & chops off some of predictions
  #bbox <- raster::extent(-180,180,-40,75)
  bbox <- raster::extent(-180,180,-50,80)

  cropnum <- 0
  for( crop in cropnames)
  {
    cropnum <- cropnum +1

    # the mirca data has zeroes in the ocean
    # I can replace 0 with NA for plotting
    # this might influence the false negative testing

    # mirca
    rst <- raster::subset(st,paste0(crop,'_mirca'))
    # for 0/1 plots
    rst[rst >0] <- 1 
    rst[rst==0] <- NA

    par(mar=c(0,0,2,0)) #bltr
    plot(rst, main=paste0(crop,' Mirca 0/1'), ext=bbox, legend=FALSE)
    plot(ne_countries(), add=TRUE, border='grey', lwd=0.1)

    # simple suit
    rst <- raster::subset(st,paste0(crop,'_suitsimp')) 
    # for 0/1 plots
    rst[rst >0] <- 1    
    rst[rst==0] <- NA

    par(mar=c(0,0,2,0)) #bltr
    plot(rst, main=paste0(crop,' simple suit'), ext=bbox, legend=FALSE)
    plot(ne_countries(), add=TRUE, border='grey', lwd=0.1)    

    # ecocrop
    rst <- raster::subset(st,paste0(crop,'_ecocrop'))
    # for 0/1 plots
    rst[rst >0] <- 1    
    rst[rst==0] <- NA

    par(mar=c(0,0,2,0)) #bltr
    plot(rst, main=paste0(crop,' ecocrop 0/1'), ext=bbox, legend=FALSE)
    plot(ne_countries(), add=TRUE, border='grey', lwd=0.1)    

  }

The new simpler method is quantitatively better than ecocrop when compared to Mirca, but still has large differences from Mirca.

#fig.height reduced to try to get 3 plots on a page

  #beware underscore in soya_bean
  cropnames <- c('maize','rice','wheat','soya_bean')

  #this doesn't currently help with plotting maps & chops off some of predictions
  #bbox <- raster::extent(-180,180,-40,75)
  bbox <- raster::extent(-180,180,-50,80)

  cropnum <- 0
  for( crop in cropnames)
  {
    cropnum <- cropnum +1

    # the mirca data has zeroes in the ocean
    # I can replace 0 with NA for plotting
    # this might influence the false negative testing

    # mirca
    rst_mirca <- raster::subset(st,paste0(crop,'_mirca'))
    # for 0/1 plots but beware of effect on calculations
    rst_mirca[rst_mirca >0] <- 1 

    # simple suit
    rst_base <- raster::subset(st,paste0(crop,'_suitsimp')) 

    rst_change <- rst_base - rst_mirca

    tm <- tm_shape(rst_change) +
      tm_raster(style='cat',
                palette = c('red4','grey95','green4'),
                labels = c('mirca only','same','baseline only'),
                title = 'suitability change') +
      tm_shape(ne_countries()) +
        tm_borders(col='grey') +
      tm_layout(legend.position = c("left","bottom"), main.title=paste0(crop,' baseline - mirca'), main.title.size=1)

    print(tm)  



  }

The utility of having simple 0/1 predictions is that it is also clearer to compare predicted change against a baseline. Calculating future - baseline gives you +1 for areas expected to gain the ability to support a crop and -1 for areas expected to lose the ability to support a crop. These can be displayed as in the map below.

maize_change <- maize_future - subset(st,"maize_suitsimp")
#plot(maize_change)


library(tmap)
tm_shape(maize_change) +
  tm_raster(style='cat',
            palette = c('red4','grey95','green4'),
            labels = c('loss','same','gain'),
            title = 'suitability change') +
  tm_shape(ne_countries()) +
    tm_borders(col='black') +
  tm_layout(legend.position = c("left","bottom"))

comparing using maximum rather than mean temp climate data

to assess simple suitability

sorghum_tmax <- ecocrop_a_raster('broom-corn',st_clim, simpler=TRUE, use_tmax=TRUE)

#test this against both tavg and Mirca
crop <- 'sorghum'
sorghum_tavg <- raster::subset(st,paste0(crop,'_suitsimp')) 
sorghum_mir <- raster::subset(st,paste0(crop,'_mirca'))

# this showed no difference at current temperatures
plot(sorghum_tmax-sorghum_tavg)
summary(sorghum_tmax-sorghum_tavg)

# but sorghum max is 40 degrees and min growth cycle is 3 months
# as Jeff said rainfall probably already excludes any areas that have >40 degrees for 3 months

# potato max is 30 degrees, does it have a difference ?
potato_tmax <- ecocrop_a_raster('potato',st_clim, simpler=TRUE, use_tmax=TRUE)
crop <- 'potato'
potato_tavg <- ecocrop_a_raster('potato',st_clim, simpler=TRUE, use_tmax=FALSE) 
#potato_mir <- raster::subset(st,paste0(crop,'_mirca'))

#aha good !
plot(potato_tavg-potato_tmax)

#compare to Mirca
folder <- 'C:\\Dropbox\\ueaHelix2017\\MIRCA\\_half_degree\\'
potato_mirca <- raster(SDMTools::read.asc.gz(paste0(folder,'annual_area_harvested_rfc_crop10_ha_30mn.asc.gz')))

#Jeff said tmax looking better, even though it does decrease the truepos rate
#but said for a paper I would then need to defend not using 
#tmin rather than tavg for min temps

#BUT is tmax better than base ecocrop for our 4 crops ?

true_pos(potato_mirca,potato_tmax)
#[1] 0.6304388
true_pos(potato_mirca,potato_tavg)
#[1] 0.7948575
true_neg(potato_mirca,potato_tmax)
#[1] 0.1350915
true_neg(potato_mirca,potato_tavg)
#[1] 0.1129034
plot(potato_mirca>0)
plot((potato_mirca>0)-potato_tmax)
plot(potato_tmax)
plot(potato_mirca>0))

# eval=FALSE for now
# Plotting maps with gradations of predictions and Mirca
# fig.height reduced to try to get 3 plots on a page

  #beware underscore in soya_bean
  cropnames <- c('maize','rice','wheat','soya_bean')

  #this doesn't currently help with plotting maps & chops off some of predictions
  #bbox <- raster::extent(-180,180,-40,75)
  bbox <- raster::extent(-180,180,-50,80)

  cropnum <- 0
  for( crop in cropnames)
  {
    cropnum <- cropnum +1

    # the mirca data has zeroes in the ocean
    # I can replace 0 with NA for plotting
    # this might influence the false negative testing

    # mirca
    rst <- raster::subset(st,paste0(crop,'_mirca'))
    rst[rst==0] <- NA

    par(mar=c(0,0,2,0)) #bltr
    plot(rst, main=paste0(crop,' Mirca'), ext=bbox)
    plot(ne_countries(), add=TRUE, border='grey', lwd=0.1)

    # simple suit
    rst <- raster::subset(st,paste0(crop,'_suitsimp'))    

    par(mar=c(0,0,2,0)) #bltr
    plot(rst, main=paste0(crop,' simple suit'), ext=bbox)
    plot(ne_countries(), add=TRUE, border='grey', lwd=0.1)    

    # ecocrop
    rst <- raster::subset(st,paste0(crop,'_ecocrop'))

    par(mar=c(0,0,2,0)) #bltr
    plot(rst, main=paste0(crop,' ecocrop'), ext=bbox)
    plot(ne_countries(), add=TRUE, border='grey', lwd=0.1)    

  }

length of crop growth cycles from ecocrop

# plot crop cycle min and max 

cropnames <- c('maize','wheat','rice','potato','soya bean','sugarcane')
cr_plot_cycle(cropnames)

References

temporary citations to start ref list [@Ramirez-Villegas2013] sorghum paper [@Anderson2015] differences in global crop models and maps [@Hijmans2012] diva-gis manual [@Estes2013] comparing mechanistic and empirical crop models [@Holzkamper2015a] climate change impacts on maize based on 3 different modelling approaches [@Iizumi2015a] how do weather influence croppping area and intensity [@Leff2004a] geographic distribution of major crops [@Muller2017a] global gridded crop model evaluation [@Portmann2010a] MIRCA crop maps

AndySouth/climcropr documentation built on May 20, 2019, 5:08 p.m.