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(dismo) #for ecocrop

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

Comparing ecocrop predictions with MIRCA crop maps at 0.5 degree resolution

#BEWARE this doesn't run from the Rmd 
#because of problems with the file_path to save outputs
#but this can just be temporary
#ideally I'll save results in package for use later
#so I could put this code into a script in 

# creating routines systematically to compare ecocrop & mirca for a range of crops

df_crops <- read_csv("name_eco, name_mir
maize, Maize
Sorghum (high altitude), Sorghum
potato, Potatoes")

# extent object for use in plots
ext <- raster::extent(-180,180,-40,75)

for( crop_num in 1:length(df_crops))
{

  name_eco <- df_crops$name_eco[crop_num]
  name_mir <- df_crops$name_mir[crop_num]  

  # repeat for rainfed and irrigated
  for( name_rain in c("rain", "irig") )
  {

    # run ecocrop
    # or extract pre-saved outfile
    # ec_out_potato_rain.grd
    name_out_eco <- paste0("ec_out_",name_eco,"_", name_rain,".grd" )

    #temporary have to do this way to be able to save the file to the package
    #file_path <- file.path("inst\\extdata", name_out_eco)

    #this is the way it should run from the package
    file_path <- system.file("extdata/", name_out_eco, package = "climcropr")

    rainfed <- ifelse(name_rain=="rain",TRUE,FALSE)

    # checking if outfile exists already
    if (! file.exists(file_path))
    {
      #message("can't find:",file_path)

      #if it doesn't exist, run ecocrop
      rst_eco <- ecocrop_a_raster( name_eco,
                                   st_clim_all = st_clim,
                                   rainfed = rainfed,
                                   filename = stringr::str_sub(file_path,1,-5), #remove .grd
                                   overwrite = TRUE)
    } else
    {
      #if it does exist load it
      rst_eco <- raster(file_path)
    }

    # allowing mirca to be switched on & off
    mirca <- TRUE

    if (mirca)
    {
      # load MIRCA file 
      rst_mir <- get_mirca(name_mir, rainfed = rainfed, plot=FALSE)
      # replace 0 with NA
      rst_mir[rst_mir==0] <- NA    

      # plot ecocrop & mirca maps one after other
      # ecocrop
      par(mar=c(0,0,2,0)) #bltr
      plot(rst_eco, main=paste0(name_eco," ecocrop prediction, ",name_rain), ext=ext)
      plot(ne_countries(), add=TRUE, border='grey', lwd=0.1)

      # mirca
      par(mar=c(0,0,2,0)) #bltr
      plot(rst_mir, main=paste0(name_mir," Mirca annual harvested area, ",name_rain), ext=ext)
      plot(ne_countries(), add=TRUE, border='grey', lwd=0.1)
    }

  } #end rain
} #end crop

#eval=FALSE

#levelplot from rsaterVis adds nice lat long histograms when applied to a single layer 
#but the breaks & colour palette didn't make  clear
#library(rasterVis)
#levelplot(rst_pot_mir_rain)
# seems like there are lots of zeroes in this mirca data 
# can I replace 0 with NA ?
rst <- rst_pot_mir_rain
rst[rst==0] <- NA

bbox <- raster::extent(-180,180,-40,75)

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

par(mar=c(0,0,2,0)) #bltr
plot(ec_out_potato_rain, main="potato ecococrop", ext=bbox) 
plot(ne_countries(), add=TRUE, border='grey', lwd=0.1)

soil ph suitability maps

0/1

# extent object for use in plots
ext <- raster::extent(-180,180,-40,75)

test_crops <- c('potato','broom-corn','maize','barley')

for(cropnum in 1:length(test_crops))
{
  cropname <- test_crops[cropnum]

  rst <- suit_soil_ph(cropname, plot=FALSE)

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

}