HICROP in R

HICROP input variables include topographic, edaphic, and climatic variables to identify spatial niche suitability for selected crops.

The topographic and edaphic variables are: slope, soil depth, soil drainage, and soil pH.

First, make some basic functions

Model funcitons

# Call package libraries
library(rgdal)
library(RColorBrewer)
library(raster)
library(dplyr)
library(tidyverse)
library(pbapply)


# # Local variables:
# # function to prompt user to input crop name
# readstring <- function() {
#   n <- readline(prompt="Enter your crop common name: ")
#   return(toString(n))
#   }
# 
# # PROMPTS USER TO INPUT CROP NAME
# crop <-  readstring() #Input, String

#crop_db <-  read.csv("CROPS.csv", stringsAsFactors = F) #insert csv file path here

# # Use col name of variable to index crop_db and return value
# locate <- function(var) {
#   n <- eval(parse(text=(paste0("crop_db$","P_OPTMIN","[na.omit(crop_db$NAME == crop)]"))))
#   return(n)
# }

# General suitability function
suit4 <- function(x) { ifelse(x < ABSmin, 0, # Conditional statement for values less than absolute minimum
                               ifelse(x > ABSmax, 0, # Conditional statement for values greater than absolute maximum
                                      ifelse(((x >= OPTmin) & (x <= OPTmax)),100, # Conditional statement for values greater than or equal absolute maximum AND less than or equal to optimal maximum (i.e. the ideal zone)
                                             ifelse(((x <= ABSmax) & (x > OPTmax)),((1-(x-OPTmax)/(ABSmax-OPTmax))*100), # Conditional statement for values less than or equal absolute maximum AND greater than optimal maximum
                                                    ifelse(((x >= ABSmin) & (x < OPTmin)),(x-ABSmin)/(OPTmin-ABSmin)*100,NA ) )))) } # Conditional statement for values greater than or equal absolute minimum AND less than optimal maximum

import crop db

# write.csv(crop_db, "crop_database_07_28_2019.csv")
# 
#  crop_dba <- crop_db
# 
# crop_db <- read.csv("crop_database_07_28_2019.csv", stringsAsFactors = FALSE)
# 
# crop_db$X <- NULL

Import crop database

Select crop

# crop defining tool
crop_define <- function() {
  n <- readline(prompt="Enter crop name: ")
  m <- crop_db
  ind <- grep(pattern = n,x = m$species,ignore.case = TRUE)
  ind <- ifelse(length(ind) == 0, grep(pattern = n,x = m$Common_names,ignore.case = TRUE), ind)
  if(length(ind) == 1) {
    print(paste0("You have selected: ", m[ind,]$species, " AKA ", strsplit(x = m[ind,67],", ")[[1]][1]))
    n <- ind
    return(select(m[n,], crop_code,NAME, species, Common_names, W_OPT_min:C_MAX,Habit, Life.form, Life.span))
  } else if(length(ind) > 1) {
    print("There are multiple potential crops. Please reveiw the list")
    print(m[ind,] %>% select(species, Common_names))
    n <- readline(prompt="There are multiple potential crops. Please reveiw the list and enter the desired CROP CODE: ")
    print(paste0("You have selected: ", m[n,]$species, " AKA ", strsplit(x = m[n,67],", ")[[1]][1]))
    return(select(m[n,], crop_code,NAME, species, Common_names, W_OPT_min:C_MAX,Habit, Life.form, Life.span))
  } else {
    if(length(ind) == 0) print(paste0(n," is not in the crop database."))
  }
}

crop_vars <- crop_define()
## [1] "You have selected: Artocarpus altilis AKA breadfruit"
crop <- crop_vars$NAME
#View(crop_vars)

#rm(crop_vars)

Start timer

# Start the clock!
ptm <- proc.time()

Define crop parameters & output file

setwd("/Users/hh/Documents/GitHub/HICROP/Outputs/")
## Error in setwd("/Users/hh/Documents/GitHub/HICROP/Outputs/"): cannot change working directory
make.dir <- function(fp) {
if(!file.exists(fp)) {  # If the folder does not exist, create a new one
make.dir(dirname(fp))
dir.create(fp)
} else {   # If it existed, delete and replace with a new one  
unlink(fp, recursive = TRUE)
dir.create(fp)
  }
} 

make.dir(paste0(crop,sep="_",Sys.Date()))
## Warning in dir.create(fp): '.' already exists
result_path <- paste0("Outputs/",crop,sep="_",Sys.Date())

Import Raster Files

setwd("/Users/hh/Documents/GitHub/HICROP")

Save crop parameters to output file

# save crop parameters to file
write.csv(crop_vars, paste0(paste0(result_path, "/params.csv")))
## Warning in file(file, ifelse(append, "a", "w")): cannot open file 'Outputs/
## breadfruit_2019-07-29/params.csv': No such file or directory
## Error in file(file, ifelse(append, "a", "w")): cannot open the connection

Ph Suitability

PH <- raster("inputs/toMac/ph.tif")
## Error in .local(.Object, ...) :
## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. (file does not exist)
# get variables
OPTmin <- crop_db$P_OPTMIN[crop_db$crop_code == crop_vars$crop_code]
OPTmax <- crop_db$P_OPTMAX[crop_db$crop_code == crop_vars$crop_code]
ABSmin <- crop_db$P_ABSMIN[crop_db$crop_code == crop_vars$crop_code]
ABSmax <- crop_db$P_ABSMAX[crop_db$crop_code == crop_vars$crop_code]

print(paste("Optimal soil pH range for", crop,"is:", OPTmin,"-",OPTmax))
## [1] "Optimal soil pH range for breadfruit is: 5.5 - 6.5"
print(paste("Absolute soil pH range for", crop,"is:", ABSmin,"-",ABSmax))
## [1] "Absolute soil pH range for breadfruit is: 4.3 - 8.7"
# calculate suitability 
P_suit <- calc(PH, suit4)
## Error in calc(PH, suit4): object 'PH' not found
# plot output
plot(P_suit, main= paste("Ph Suitability for ",crop), breaks = c(0,25,50,75,100), col = rainbow(5),
     sub=paste("Optimal pH range:", OPTmin,"-",OPTmax, "\n Absolute pH range:", ABSmin,"-",ABSmax))

plot of chunk ph Suitability

# cleanup & save
writeRaster(P_suit, filename = paste0("Outputs/",crop,sep="_",Sys.Date(),"/ph_suit.tif"))
## Error in .local(x, filename, ...): Attempting to write a file to a path that does not exist:
##   Outputs/breadfruit_2019-07-29
rm(PH)
## Warning in rm(PH): object 'PH' not found

A primary difference between the two is that the python script performs multiple condtional operations (i.e., rastcon1 through rastcon5) and output operations (i.e., true3 and true4) for each input value, before determining which conditional and output apply (see rast_suit, line 138). The R code reduces the total number of calculations by only determining outputs once a conditional is met (see soilsuit, lines 202-206). This truncating approach is taken for all scripts below.

Other than that the only differences are the that the python scripts pull in parameter data from a csv whereas currently only dummy data was selected and hardcorded in (e.g., ph_ABSmax etc).

Depth Suitability

DEPTH <- raster("inputs/toMac/depths.tif")
## Error in .local(.Object, ...) :
## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. (file does not exist)
#freq(DEPTH)

# Get variables 
OPTmin <- crop_db$D_OPT[crop_db$crop_code == crop_vars$crop_code]
ABSmin <- crop_db$D_ABS[crop_db$crop_code == crop_vars$crop_code]

print(paste("Optimal minimum soil depth for", crop,"is:", OPTmin, "cm"))
## [1] "Optimal minimum soil depth for breadfruit is: 150 cm"
print(paste("Absolute minimum soil depth for", crop,"is:", ABSmin, "cm"))
## [1] "Absolute minimum soil depth for breadfruit is: 20 cm"
# # define custom function
# suit2min <- function(x) { ifelse(x < ABSmin, 0, # Conditional statement for values less than absolute minimum
#                                  ifelse(x >= OPTmin, 100, # Conditional for values greater than or equal to optimal minimum
#                                         ifelse(((x >= ABSmin) & (x < OPTmin)), ((x-ABSmin)/(OPTmin-ABSmin)*100), NA)))} 
# 
# # run calculation
# D_suit <- calc(DEPTH, suit2min)

# use pre-run suitability stack
D_suit <- depth_stack[[which(apply(dep, 1, function(r) any(r[1] %in% OPTmin & r[2] %in% ABSmin)))]]


# plot output
plot(D_suit, main= paste("Depth Suitability for ",crop), breaks = c(0,25,50,75,100), col = rainbow(5),
     sub = paste('Optimal minimum depth is ', OPTmin, 'cm, \n Absolute minimum depths is', ABSmin, 'cm.'))

plot of chunk Depth Suitability

# save and cleanup
writeRaster(D_suit, filename = paste0("Outputs/",crop,sep="_",Sys.Date(),"/depth_suit.tif"))
## Error in .local(x, filename, ...): Attempting to write a file to a path that does not exist:
##   Outputs/breadfruit_2019-07-29
#rm(DEPTH)

Drainage Suitability

#DRAIN <- raster("inputs/toMac/drain.tif")
#freq(DRAIN)

# get variables
OPTmin <- crop_db$W_OPT_min[crop_db$crop_code == crop_vars$crop_code]
OPTmax <- crop_db$W_OPT_max[crop_db$crop_code == crop_vars$crop_code]
ABSmin <- crop_db$W_ABS_min[crop_db$crop_code == crop_vars$crop_code]
ABSmax <- crop_db$W_ABS_max[crop_db$crop_code == crop_vars$crop_code]

print(paste("Optimal soil drainage code range for", crop,"is:", OPTmin,"-",OPTmax))
## [1] "Optimal soil drainage code range for breadfruit is: 6 - 6"
print(paste("Absolute soil drainage code range for", crop,"is:", ABSmin,"-",ABSmax))
## [1] "Absolute soil drainage code range for breadfruit is: 6 - 6"
# calculate suitability
#W_suit <- calc(DRAIN, suit4)

# use pre-run suitability stack
W_suit <- drain_stack[[which(apply(dra, 1, function(r) any(r[1] %in% OPTmin & r[2] %in% OPTmax & r[3] %in% ABSmin & r[4] %in% ABSmax)))]]

plot(W_suit, main= paste("Drainage Suitability for ",crop), breaks = c(0,25,50,75,100), col = rainbow(5), 
     sub=paste("Optimal drainage range:", OPTmin,"-",OPTmax, "\n Absolute drainage range:", ABSmin,"-",ABSmax))

plot of chunk Drainage Suitability

# save and cleanup
writeRaster(W_suit, filename = paste0("Outputs/",crop,sep="_",Sys.Date(),"/drainage_suit.tif"))
## Error in .local(x, filename, ...): Attempting to write a file to a path that does not exist:
##   Outputs/breadfruit_2019-07-29
#rm(DRAIN)

Slope Suitability

# SLOPE <- raster("inputs/toMac/slope.tif")

# Determine slope input variables
OPTmax <- crop_db$S_OPTMAX[crop_db$crop_code == crop_vars$crop_code]
ABSmax <- crop_db$S_ABSMAX[crop_db$crop_code == crop_vars$crop_code]

# Reivew texts
print(paste("Optimal maximum slope for", crop,"is:", OPTmax, "%"))
## [1] "Optimal maximum slope for breadfruit is: 28 %"
print(paste("Absolute maximum soil slope for", crop,"is:", ABSmax, "%"))
## [1] "Absolute maximum soil slope for breadfruit is: 38 %"
# 
# # Use fuzzy calculation to determine range of suitability for
# suit2max <- function(x) { ifelse(x > ABSmax, 0, # Conditional statement for values greater than absolute maximum slope
#                                   ifelse(x <= OPTmax, 100, # Conditional statement for values less than or equal to optimal maximum slope
#                                          ifelse(((x <= ABSmax) & (x > OPTmax)),((1-(x - OPTmax)/(ABSmax-OPTmax))*100), NA)))} # Conditional statement for values greater than or equal absolute minimun AND less than optimal minimum
# 
# # calculate suitability
# S_suit <- calc(SLOPE, suit2max)

#Plot your output
plot(S_suit, main= paste("Slope Suitability for ",crop), breaks = c(0,25,50,75,100), col = rainbow(5), 
     sub = paste('Optimal maximum slope is ', OPTmax, '%, \n Absolute maximum slope is', ABSmax, '%.'))

plot of chunk Slope Suitability

# save and cleanup
writeRaster(S_suit, filename = paste0("Outputs/",crop,sep="_",Sys.Date(),"/slope_suit.tif"))
## Error in .local(x, filename, ...): Attempting to write a file to a path that does not exist:
##   Outputs/breadfruit_2019-07-29
#rm(SLOPE)

Calculate Temperature

# Get temp rasters
# get file list
# read in as rasters
# store as raster stack

Env_Tmax <- stack(sapply(Sys.glob("inputs/toMac_climate/tmax/*.tif"), raster))
## Error in x[[1]]: subscript out of bounds
Env_Tmin <- stack(sapply(Sys.glob("inputs/toMac_climate/tmin/*.tif"), raster))
## Error in x[[1]]: subscript out of bounds
# Get float values for temperature variables from ecocrop db
ABSmax <- crop_db$T_ABSMAX[crop_db$crop_code == crop_vars$crop_code]
ABSmin <- crop_db$T_ABSMIN[crop_db$crop_code == crop_vars$crop_code]
OPTmax <- crop_db$T_OPTMAX[crop_db$crop_code == crop_vars$crop_code]
OPTmin <- crop_db$T_OPTMIN[crop_db$crop_code == crop_vars$crop_code]


print(paste("Optimal Temperature range for", crop,"is:", OPTmin,"-",OPTmax, "°C"))
## [1] "Optimal Temperature range for breadfruit is: 21 - 33 °C"
print(paste("Absolute Temperature range for", crop,"is:", ABSmin,"-",ABSmax, "°C"))
## [1] "Absolute Temperature range for breadfruit is: 16 - 40 °C"
## Suitability Evaluation
Tmax <- calc(Env_Tmax, suit4)
## Error in .local(.Object, ...):
Tmin <- calc(Env_Tmin, suit4)
## Error in .local(.Object, ...):
# make a raster stack with minimum value from min and max suitability layers by month

f <- function(x1, x2)  ifelse(x1 < x2, x1, x2)
T_suit <- overlay(Tmax, Tmin, fun=f, forcefun=TRUE)


#plot(T_suit, main= paste0("Temperature Suitability for ", crop), breaks = c(0,25,50,75,100), col = rainbow(5))

library(rasterVis)

levelplot(T_suit, main= paste0("Temperature Suitability for ", crop), 
          names.attr=month.name[1:12], 
          sub=paste("Optimal Temperature range:", OPTmin,"-",OPTmax, "°C \n Absolute Temperature range:", ABSmin,"-",ABSmax, "°C"),
          at = seq(0, 100, by = 20),
          col.regions = rainbow(5))

plot of chunk Temperature Suitability

Rainfall suitability analysis

# get list of Tif filenames and make raster stack
# NEED TO FIGURE OUT WHY THIS BROKE
# Env_rain <- sapply(Sys.glob("inputs/rainfall/*.tif"), raster) %>% stack() 
Env_rain <- stack(sapply(Sys.glob("inputs/toMac_climate/rain/*.tif"), raster))
## Error in x[[1]]: subscript out of bounds
# get variables
ABSmax <- round(crop_db$R_ABSMAX[crop_db$crop_code == crop_vars$crop_code]/12,2)
ABSmin <- round(crop_db$R_ABSMIN[crop_db$crop_code == crop_vars$crop_code]/12,2)
OPTmax <- round(crop_db$R_OPTMAX[crop_db$crop_code == crop_vars$crop_code]/12,2)
OPTmin <- round(crop_db$R_OPTMIN[crop_db$crop_code == crop_vars$crop_code]/12,2)

# Get float values for precipitation variables from ecocrop csv file:
print(paste("Optimal Precipitation range for", crop,"is:", OPTmin,"-",OPTmax, "mm"))
## [1] "Optimal Precipitation range for breadfruit is: 125 - 250 mm"
print(paste("Absolute Precipitation range for", crop,"is:", ABSmin,"-",ABSmax, "mm"))
## [1] "Absolute Precipitation range for breadfruit is: 83.33 - 291.67 mm"
# calculate suitability
R_suit <- calc(Env_rain, suit4)
## Error in .local(.Object, ...):
# plot outputs
levelplot(R_suit, main= paste0("Precipitation Suitability for ", crop), 
          names.attr=month.name[1:12], 
          sub=paste("Optimal Precipitation range:", round(OPTmin),"-",round(OPTmax), "mm \n Absolute Precipitation range:", round(ABSmin),"-",round(ABSmax), "mm"),
          at = seq(0, 100, by = 20),
          col.regions = rainbow(5))

plot of chunk Rainfall Suitability

Seasonal Calculations

# Get the max days of season, 
# divide and round up to produce integer for season length in months

#Cmin <- as.integer(ceiling(locate("C_MIN")/30.42)) 
Cmax <- as.integer(ceiling(crop_db$C_MAX[crop_db$crop_code == crop_vars$crop_code]/30.42)) # 365/12=30.42

print(paste("Maximum growing season length for", crop,"is:", Cmax, "months"))
## [1] "Maximum growing season length for breadfruit is: 12 months"
#MAKE ONE GIANT SEASON ASSESSMENT FUNCTION

# Calculate Average Suitability and Difference across months in a season
season_avg <- function(var){ 
  # make a double stack of monthly input variable (R or T)
  # this will allow for calculating across all season lenghts (e.g., 12 month season starting in December..)
  months24 <- stack(var, var)

  # Make a 2-row matrix of all growing season possibiliites, based on maximum season length
  # index season by col
  seasons <- sapply(1:12, function(n) c(n, n+Cmax-1))
  # make a list of stack 
  z <- sapply(1:12, function(i) months24[[seasons[1,i]:seasons[2,i]]])

# calculate seasonal means
agg.fun <- function(x,...) { ## you should add ... for na.rm,and other extra parameters
   # if Cmax >= 12, run only the first season stack mean, else run mean for each season
  if (Cmax >= 12){return(calc(x[1][[1]], mean)) # if crop cycle lengths is 12 or more months, take the mean of year suitability
  print("Calculating mean season")
  }else{ # else, make mean of all seasons
    print("Calculating mean seasons")
    return(stack(pbsapply(1:12, function(i) calc(z[i][[1]], mean))))  # REPLACE wih sapply(1:12,... and remove stack?
  }
  } 

 seasons_mean <- agg.fun(z)
  return(seasons_mean)
}


#Check if RAINFALL / TEMP is annual or monthly
Rcat <- if (dim(R_suit)[3] == 12){
  print("Evaluating RAINFALL suitability by season")
  season_avg(R_suit) 
}
## [1] "Evaluating RAINFALL suitability by season"
Tcat <- if (dim(T_suit)[3] == 12){
  print("Evaluating TEMPERATURE suitability by season")
  season_avg(T_suit) 
}
## [1] "Evaluating TEMPERATURE suitability by season"
#Create Rain by Temp suitability 
rainxtemp <- function(x,y){
  rxt_stack <- stack() # make and empty raster stack
  rast_stack <- stack(x,y) # stack two input raster stacks (Rcat & Tcat) as one
 if(dim(rast_stack)[3] == 24){ # if season is less than 12 months
        print("Calculating Rain by Temp suitability across multiple seasons")
   for(i in 1:12){ # for each season mean pairing
          print(paste0("Calculating Rain by Temp suitability for season ", i))
     rxt <- (rast_stack[[i]] * rast_stack[[i + 12]]) / 100 # factor season mean pairs (Rcat*Tcat)
          rxt_stack <- stack(rxt_stack, rxt)} # store result in new stack
          }
  else { #if Cmax >= 12, all seasons will be the same (i.e. there is no need to loop through additional seasons)
    print("Calculating Rain by Temp for complete year")
    fun <- function(z) { ((z[[1]] * z[[2]]) / 100)}
    rxt_stack <- calc(rast_stack, fun) # store result in new stack
    }
  return(rxt_stack) # return new stack
}


#if TEMP & RAIN is monthly, then run the rainxtemp function defined above
RxxT <- if((dim(R_suit)[3] == 12) & (dim(T_suit)[3] == 12)) rainxtemp(Rcat, Tcat)
## [1] "Calculating Rain by Temp for complete year"
# plot and save if Cmax >=12 aka 12+ month season
if (dim(RxxT)[3]==1) plot(RxxT)

plot of chunk Seasonal Calculations

if (dim(RxxT)[3]==1) writeRaster(RxxT, paste0(result_path, "/RxT_suit.tif"))
## Error in .local(x, filename, ...): Attempting to write a file to a path that does not exist:
##   Outputs/breadfruit_2019-07-29

Final Suitability Overlay & 0ptimal Planting Month

if (dim(T_suit)[3] == 12) print("Temp is monthly (seasonal")
## [1] "Temp is monthly (seasonal"
if (dim(T_suit)[3] == 1) print("Temp is annual")
if (dim(R_suit)[3] == 12) print("Precip is monthly (seasonal")
## [1] "Precip is monthly (seasonal"
if (dim(R_suit)[3] == 1) print("Precip is annual")




#if either rainfall or temperature is seasonal {ANNUAL?}, determine optimal season
if (any((dim(R_suit)[3] == 12),(dim(T_suit)[3] == 12))) {
  #if both rainfall and temperature are monthly 
  # use RxT raster catalog created by "Seasons" script
  if ((dim(R_suit)[3] == 12) & (dim(T_suit)[3] == 12)) {
    seasonal <- RxxT
  } else if((dim(R_suit)[3] == 1) & (dim(T_suit)[3] == 12)) {
    seasonal <- T_suit
  } else if((dim(R_suit)[3] == 12) & (dim(T_suit)[3] == 1)) {
    seasonal <- R_suit
  } 

#Determine Optimal Season
# make function to find index of max cell in stack
which.max.na <- function(x, ...) ifelse(length(x) == sum(is.na(x)), 0, which.max(x))
# run
opt_season <- calc(seasonal, which.max.na)
#set all pixels that are not contained in month_nums to NA
opt_season[!(opt_season %in% c(1:12))] <- NA



#Return MAXIMUM Score (from optimal season)
  seasuit <- max(seasonal)
  print ("Maximum Seasonal (Temp and/or Rain) Suitability Score Raster: seasuit")

} else print("All criteria are annual")
## Warning in max(new("RasterLayer", file = new(".RasterFile", name = "/
## private/var/folders/xw/8jnc9qm53jz_gn09wx6tc3yr0000gn/T/RtmpWAr4ia/raster/
## r_tmp_2019-07-29_113539_42736_45557.grd", : Nothing to summarize if you
## provide a single RasterLayer; see cellStats
## [1] "Maximum Seasonal (Temp and/or Rain) Suitability Score Raster: seasuit"
# # RESAMPLING TO BE REMOVED
# Psuit <- resample(P_suit,seasuit,method='ngb')
# Dsuit <- resample(D_suit,seasuit,method='ngb')
# Ssuit <- resample(S_suit,seasuit,method='ngb')
# Wsuit <- resample(W_suit,seasuit,method='ngb')
# 
# 
# #Overall MAXIMUM Score (all criteria included, across all available seasons)
# if ((dim(R_suit)[3] == 12) & (dim(T_suit)[3] == 12)) {
#   all_criteria <- stack(Dsuit, Psuit, Ssuit, Wsuit, seasuit)
# } else if((dim(R_suit)[3] == 1) & (dim(T_suit)[3] == 12)) {
#   all_criteria <- stack(Dsuit, Psuit, Ssuit, Wsuit, seasuit, annual_rainfall)
# } else if((dim(R_suit)[3] == 12) & (dim(T_suit)[3] == 1)) {
#   all_criteria <- stack(Dsuit, Psuit, Ssuit, Wsuit, seasuit, annual_temp)
# } else if ((dim(R_suit)[3] == 1) & (dim(T_suit)[3] == 1)) {
#   all_criteria <- stack(Dsuit, Psuit, Ssuit, Wsuit, annual_rainfall, annual_temp)
# } else print("error in at least one raster catalog, catalog should contain 1 or 12 rasters")
# 


#Overall MAXIMUM Score (all criteria included, across all available seasons)
if ((dim(R_suit)[3] == 12) & (dim(T_suit)[3] == 12)) {
  all_criteria <- stack(D_suit, P_suit, S_suit, W_suit, seasuit)
} else if((dim(R_suit)[3] == 1) & (dim(T_suit)[3] == 12)) {
  all_criteria <- stack(D_suit, P_suit, S_suit, W_suit, seasuit, annual_rainfall)
} else if((dim(R_suit)[3] == 12) & (dim(T_suit)[3] == 1)) {
  all_criteria <- stack(D_suit, P_suit, S_suit, W_suit, seasuit, annual_temp)
} else if ((dim(R_suit)[3] == 1) & (dim(T_suit)[3] == 1)) {
  all_criteria <- stack(D_suit, P_suit, S_suit, W_suit, annual_rainfall, annual_temp)
} else print("error in at least one raster catalog, catalog should contain 1 or 12 rasters")



suit <- min(all_criteria)

plot(suit, main= paste("Suitability for ",crop), breaks = c(0,25,50,75,100), col = rainbow(5))

plot of chunk Final Suitability

# cleanup & save
writeRaster(suit, filename = paste0("Outputs/",crop,sep="_",Sys.Date(),"/all_suit.tif"))
## Error in .local(x, filename, ...): Attempting to write a file to a path that does not exist:
##   Outputs/breadfruit_2019-07-29
writeRaster(opt_season, filename = paste0("Outputs/",crop,sep="_",Sys.Date(),"/optimal_season.tif"))
## Error in .local(x, filename, ...): Attempting to write a file to a path that does not exist:
##   Outputs/breadfruit_2019-07-29

Limiting Factors

print(paste("Maximum growing season length for", crop,"is:", Cmax, "months"))
## [1] "Maximum growing season length for breadfruit is: 12 months"
# ######
# #use a seasonal raster if it exists otherwise use a non-seasonal raster
# Rsuit_cat <- env.workspace +"/"+ crop[0:3] + "_Rsuit" + str(Cmax)
# if not arcpy.Exists(Rsuit_cat):
#   R_Rsuit = env.workspace +"/"+ crop[0:3] + "_Rsuit/Raster.OBJECTID=1"
# 
# Tsuit_cat = env.workspace +"/"+ crop[0:3] + "_Tsuit" + str(Cmax)
# if not arcpy.Exists(Tsuit_cat):
#   R_Tsuit = env.workspace +"/"+ crop[0:3] + "_Tsuit/Raster.OBJECTID=1"
# ################

# If season is >= 12 months uses mean of R_suit stack (Rcat), else uses mean of first season stack (Rcat[[1]])
### CHECK IF THESE TWO LINES ARE NEEDED
#Rc <- ifelse (dim(Rcat)[3] == 12, Rcat[[1]], Rcat)
#Rc <- ifelse (dim(Rcat)[3] == 1, Rcat, Rcat[[1]])


if (dim(Rcat)[3] == 12) {
  Rc <- Rcat[[1]]
} else {
  Rc <-Rcat
}
Rc
## class      : RasterLayer 
## dimensions : 3663, 5712, 20923056  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 370100, 941300, 2093300, 2459600  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=4 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : /private/var/folders/xw/8jnc9qm53jz_gn09wx6tc3yr0000gn/T/RtmpWAr4ia/raster/r_tmp_2019-07-29_113203_42736_43794.grd 
## names      : layer 
## values     : 0, 90.93108  (min, max)
if (dim(Tcat)[3] == 12) {
  Tc <- Tcat[[1]]
} else {
  Tc <-Tcat
}
Tc
## class      : RasterLayer 
## dimensions : 3663, 5712, 20923056  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 370100, 941300, 2093300, 2459600  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=4 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : /private/var/folders/xw/8jnc9qm53jz_gn09wx6tc3yr0000gn/T/RtmpWAr4ia/raster/r_tmp_2019-07-29_113354_42736_32893.grd 
## names      : layer 
## values     : 0, 99.96182  (min, max)
ss <- stack(P_suit, D_suit, W_suit, S_suit, Tc, Rc) # make stack of layers
names(ss) <- c("Ph","Depth","Drainage","Slope","Temp","Rain") # rename layers for no reason
mns <- whiches.min(ss) # get lowest value(s) by layer

#freq(mns) %>% View()
#plot(mns)
#####

#unique(mns)

# access raster attribute table
mns <- ratify(mns) # make a rat
#mns
rat <- levels(mns)[[1]] 
#rat

# make new column and substitue layer index num (1:6) for layer name strings
rat$limit <- gsub(1, "PH_", rat$ID)
rat$limit <- gsub(2, "DEPTH_", rat$limit)
rat$limit <- gsub(3, "DRAIN_", rat$limit)
rat$limit <- gsub(4, "SLOPE_", rat$limit)
rat$limit <- gsub(5, "TEMP_", rat$limit)
rat$limit <- gsub(6, "RAIN_", rat$limit)
rat
##        ID                           limit
## 1       1                             PH_
## 2       2                          DEPTH_
## 3       3                          DRAIN_
## 4       4                          SLOPE_
## 5       5                           TEMP_
## 6       6                           RAIN_
## 7      12                       PH_DEPTH_
## 8      13                       PH_DRAIN_
## 9      14                       PH_SLOPE_
## 10     15                        PH_TEMP_
## 11     16                        PH_RAIN_
## 12     23                    DEPTH_DRAIN_
## 13     24                    DEPTH_SLOPE_
## 14     25                     DEPTH_TEMP_
## 15     26                     DEPTH_RAIN_
## 16     34                    DRAIN_SLOPE_
## 17     35                     DRAIN_TEMP_
## 18     36                     DRAIN_RAIN_
## 19     45                     SLOPE_TEMP_
## 20     46                     SLOPE_RAIN_
## 21     56                      TEMP_RAIN_
## 22    123                 PH_DEPTH_DRAIN_
## 23    124                 PH_DEPTH_SLOPE_
## 24    125                  PH_DEPTH_TEMP_
## 25    126                  PH_DEPTH_RAIN_
## 26    134                 PH_DRAIN_SLOPE_
## 27    135                  PH_DRAIN_TEMP_
## 28    136                  PH_DRAIN_RAIN_
## 29    145                  PH_SLOPE_TEMP_
## 30    146                  PH_SLOPE_RAIN_
## 31    156                   PH_TEMP_RAIN_
## 32    234              DEPTH_DRAIN_SLOPE_
## 33    235               DEPTH_DRAIN_TEMP_
## 34    236               DEPTH_DRAIN_RAIN_
## 35    245               DEPTH_SLOPE_TEMP_
## 36    246               DEPTH_SLOPE_RAIN_
## 37    256                DEPTH_TEMP_RAIN_
## 38    345               DRAIN_SLOPE_TEMP_
## 39    346               DRAIN_SLOPE_RAIN_
## 40    356                DRAIN_TEMP_RAIN_
## 41    456                SLOPE_TEMP_RAIN_
## 42   1234           PH_DEPTH_DRAIN_SLOPE_
## 43   1235            PH_DEPTH_DRAIN_TEMP_
## 44   1236            PH_DEPTH_DRAIN_RAIN_
## 45   1245            PH_DEPTH_SLOPE_TEMP_
## 46   1246            PH_DEPTH_SLOPE_RAIN_
## 47   1256             PH_DEPTH_TEMP_RAIN_
## 48   1345            PH_DRAIN_SLOPE_TEMP_
## 49   1346            PH_DRAIN_SLOPE_RAIN_
## 50   1356             PH_DRAIN_TEMP_RAIN_
## 51   1456             PH_SLOPE_TEMP_RAIN_
## 52   2345         DEPTH_DRAIN_SLOPE_TEMP_
## 53   2346         DEPTH_DRAIN_SLOPE_RAIN_
## 54   2356          DEPTH_DRAIN_TEMP_RAIN_
## 55   2456          DEPTH_SLOPE_TEMP_RAIN_
## 56   3456          DRAIN_SLOPE_TEMP_RAIN_
## 57  12345      PH_DEPTH_DRAIN_SLOPE_TEMP_
## 58  12346      PH_DEPTH_DRAIN_SLOPE_RAIN_
## 59  12356       PH_DEPTH_DRAIN_TEMP_RAIN_
## 60  13456       PH_DRAIN_SLOPE_TEMP_RAIN_
## 61  23456    DEPTH_DRAIN_SLOPE_TEMP_RAIN_
## 62 123456 PH_DEPTH_DRAIN_SLOPE_TEMP_RAIN_
levels(mns) <- rat # assign levels based on new rat

#freq([email protected]@attributes$limit)

#mns

writeRaster(mns,paste0(result_path,"/limits.tif"))
## Error in .local(x, filename, ...): Attempting to write a file to a path that does not exist:
##   Outputs/breadfruit_2019-07-29

Validation: Ommision Rate & RMSE

# Name:         RMSE
# Purpose:      Calculate root mean square error based on known crop locations

library(dplyr)
library(raster)
# make shapefile of 
baseline <- shapefile("~/GIS/GIS Files/2015AgBaseline.shp/2015AgBaseline.shp")

bs_crops <- [email protected]$CropCatego %>% unique()

# check for crop in location file, save as raster if present
if (crop %in% bs_crops) { 
  print(paste0(crop," has location data. RMSE & OR will be calculated."))
  crop_loc <- baseline %>% subset(., CropCatego %in% crop) %>% #subset shapefile to crop data
    rasterize(., seasuit) #rasterize to other layer parameters
  crop_loc[crop_loc >= 1] <- 1 # set values to 1
  suitscore <- suit
  # Calculate Root Mean Square Error (RMSE) for everywhere
  rmse <- ((suitscore/100)-1)**2
  # trim data to only areas with crop location
  rmse <- mask(rmse, crop_loc)
  # get sum of RMSE values
  rmse_sum <- cellStats(x = rmse,stat = 'sum')
  # count cells with crop 
  rmse_count <- zonal(rmse, crop_loc, 'count')
  # calc RMSE value
  rmse_value <- rmse_sum/rmse_count[1,2]
  # report
  print(paste0(crop," RMSE = ", round(rmse_value,5)))
  # Calculate Omission Rate (OR)
  #reclassify values below above 0 to 1
  suitscore[suitscore>0]=1
  suitscore <- mask(suitscore, crop_loc)  # trim to just location data
  or_sum <- cellStats(x = suitscore,stat = 'sum') # get sum of viable cells with crop present
  or_value =  1-(or_sum/rmse_count[1,2])
  print(paste0(crop," Omission Rate= ", round(or_value,5)))  # report
  # create a table
  rmse_table <- data.frame("Crop" = crop, "Count" = rmse_count[1,2] , "SUM_RMSE" = rmse_sum, "RMSE" = rmse_value, "SUM_OR" = or_sum, "OMISSION" = or_value)
print(rmse_table)
#save outputs & cleanup
print("Saving RMSE & OR outputs")
writeRaster(rmse, filename = paste0("Outputs/",crop,sep="_",Sys.Date(),"/rmse.tif"))
write.csv(rmse_table, paste0("Outputs/",crop,sep="_",Sys.Date(),"/RMSE_OR.csv"))
rm(suitscore, rmse, rmse_count, rmse_value, or_sum, or_value)
  }else  {print(paste0("No location data for ",crop, ". RMSE & Omission Rate calculated for HDOA 2015 Baseline Study crops only"))
}
## [1] "No location data for breadfruit. RMSE & Omission Rate calculated for HDOA 2015 Baseline Study crops only"
# PAU!

Stop the clock and Tell me its over

# Stop the clock
print(paste0(round((proc.time() - ptm)[3]/60,0), " minute runtime."))
## [1] "9 minute runtime."
notify(paste0("R is done, after a ", round((proc.time() - ptm)[3]/60,0), " minute runtime."), speakIt=TRUE)
## Warning in system(cmd): error in running command