Grass-Cast in the Southwest

An analysis of its forecast skill

Introduction

Grass, forbs, and Shrubs. Rangelands cover more than 54% of the world’s surface, performing key roles in carbon storage, and providing 70% of the world’s livestock feed source. (Brown & Thorpe, 2008). Rangelands are among the most threatened and forsaken ecosystems in the USA, with less than 10% of climate goals stirring towards their protection and over 2.1 million acres lost just in 2019 (ILRI, 2021). Among rangelands, the American Southwest has unique challenges of augmenting wildfire occurrence and increasing aridity due to climate change (Seager et al., 2007). Efficacy in the management, among other conservation tactics should be used to preserve the desert rangelands of the southwest USA. Among these efforts is Grass-Cast. This tool can predict rangelands’ biomass using 30 years of satellite data to prevent overgrazing and plan actions for subsequent seasons. Assessment of the accuracy of the prediction is focused on the Great plains in the USA (Hartman et al., 2020), but further information for the efficacy in the southwest is still needed.

The following project uses data available in Grass-Cast Explorer. I will compare weekly data of predicted ANPP (above-ground net primary productivity) and NDVI (Normalized difference vegetation index) from each forecast during the growing season and compare it to the observable ANPP. Then, we will choose an appropriate statistical model to assess the predictability of the tool, with satellite data and on-ground data. We hope that this project will support the team’s work behind Grass Cast by giving them a fresh perspective on the tool forecast capabilities for the southwest. We will also hope that ranchers and government agencies find in Grass Cast a mature tool that will help them plan the use and conservation of rangelands across the southwest.

Forecast

The forecast analysis of Grass-Cast is based on the paper by Hartman et al, and is strongly correlated with the the weather forecast and history of NOAA(National Oceanic and Atmospheric Administration), mostly the precipitation (Figure 1). This input is add in to the DayCent model, and ecosystem-level biochemical model that simulates various ecological process, the main value extracted from this model is evapotranspiration (AET). Later, the model takes information of NDVI (Normalized Difference Vegetation Index) from the MODIS sensor in the Terra satellite, NDVI measures the greenness or the the state of the vegetation, where higher levels of NDVI means healthier vegetation. Finally, three ANPP (Above ground net primary production) scenarios are computed (average, above average and bellow average) measuring the biomass in pounds per acre.

Figure1

Data

The forecast data is from the Grass-Cast platform can be accessed online (National Drought Mitigation Center, 2022).

Methodology

For the southwest, there are two growing seasons Spring (April-May) and Summer(August-September), the spring is mostly correlated with the winter rains (November-March) and the summer is correlated with the monsoon season (June-September), (Smith et al., 2019).

In this experiment we are trying to understand how well the forecast correlates to the final observation of Grass Cast and also see the effects of rainfall in the final forecast. To accomplish this objective I will take coefficient of determination (R2) between the percent difference in ANPP from each forecast during the growing season compared to the percent difference in ANPP for the final forecast. The result will be a line graft of a correlation from 8 forecast of 2020 and the 12 forecast of 2021 in the three different scenarios to the final observation of the year.

Beside this, I will interpolate rainfall data from NOAA from each growing season to explain the possible reasons for the data to change correlations during time.

Coding

##Charge the libraries 
#Data Manipulation
library(tidyverse) #pipes and more
library(janitor) #cleans col numbers and more
library(data.table) #for splits

#Graphics Tools
library(ggplot2) #draws the graphs
library(cowplot) #plots side by side
library(ggpubr) #adds ggarange
library(tidytext) # for the point 
library(ggrepel)
library(patchwork) #adds ggplots

#Spatial Tools
library(maps) #database of USA maps
library(sp)
library(sf)
library(raster) 
library(rgeos)

###Nearest neighbour
library(dismo) 
library(gstat)

Correlation Graph

####################Data base prep
#upload the data and select collumns 
GC <- read.csv(file = "data/az_nm_2000_2020.csv", head = TRUE, sep=",") 


gfg_data <- read.csv(file = "data/grass_cast_20-21.csv", head = TRUE, sep=",")   %>%   subset(select=c("gridID","Year","Forecast","NPP_predict_below",
"NPP_predict_avg","NPP_predict_above","Order"))

summerwinter <- read.csv(file = "data/RatioSummerWinter.csv", head = TRUE, sep=",") %>% subset(select=c("gridID","latitude","longitude","pptRatioSummerWinter")) 

#Joins the table above
RatioSW_ForecastANPP <- left_join(gfg_data, summerwinter, by = c("gridID" = "gridID"))

#split each data in summer, transision, winter
ANPP_FORECAST_W <- subset(RatioSW_ForecastANPP, pptRatioSummerWinter < '0.8')
ANPP_FORECAST_T <- RatioSW_ForecastANPP  %>%   filter(pptRatioSummerWinter >0.8 &  pptRatioSummerWinter < 1.2)
ANPP_FORECAST_S <- subset(RatioSW_ForecastANPP, pptRatioSummerWinter > '1.2')

#filter year 2020
ANPP_FORECAST_2020_W <- subset(ANPP_FORECAST_W, Year == '2020') 
ANPP_FORECAST_2020_T <- subset(ANPP_FORECAST_T, Year == '2020')
ANPP_FORECAST_2020_S <- subset(ANPP_FORECAST_S, Year == '2020')

#Slip the forecast
ANPP_FORECAST_2020_W <- split(ANPP_FORECAST_2020_W,ANPP_FORECAST_2020_W$Forecast)
ANPP_FORECAST_2020_T <- split(ANPP_FORECAST_2020_T,ANPP_FORECAST_2020_T$Forecast)
ANPP_FORECAST_2020_S <- split(ANPP_FORECAST_2020_S,ANPP_FORECAST_2020_S$Forecast)

#filter year 2021
ANPP_FORECAST_2021_W <- subset(ANPP_FORECAST_W, Year == '2021')
ANPP_FORECAST_2021_T <- subset(ANPP_FORECAST_T, Year == '2021')
ANPP_FORECAST_2021_S <- subset(ANPP_FORECAST_S, Year == '2021')

#Slip the forecast
ANPP_FORECAST_2021_W <- split(ANPP_FORECAST_2021_W,ANPP_FORECAST_2021_W$Forecast)
ANPP_FORECAST_2021_T <- split(ANPP_FORECAST_2021_T,ANPP_FORECAST_2021_T$Forecast)
ANPP_FORECAST_2021_S <- split(ANPP_FORECAST_2021_S,ANPP_FORECAST_2021_S$Forecast)

#The data set is ready 

####################################Functions and plots

#this functions takes the last date and compares to the every other forecast  and see how well is the forecast predicting 

fc <- function(n, x, y) {
  
  #N is the observation we are testing
  #x is the final observation
  #y is data frame
  
  a <- cor(y[[n]]$NPP_predict_above, y[[x]]$NPP_predict_above)
  b <- cor(y[[n]]$NPP_predict_avg, y[[x]]$NPP_predict_avg)
  c <- cor(y[[n]]$NPP_predict_below, y[[x]]$NPP_predict_below)

  #Creates a matrtix 
  ANPP.r <- matrix(c(a,b,c),ncol = 3,byrow = TRUE)
  colnames(ANPP.r) <- c('Abv','Avg','Blw')
  ANPP.r <- as.data.frame(ANPP.r) 
  #chages row names to date
  row.names(ANPP.r) <- as.Date(names(y[x]), format = "%m/%d/%Y")  
  
  return(ANPP.r)
}

##Creates an empty list for the for loop
r.list2020_W <- list()
r.list2020_T <- list()
r.list2020_S <- list()

r.list2021_W <- list()
r.list2021_T <- list()
r.list2021_S <- list()

##For loop to have all the rsquer in the list
#the final observation was on place 8 for 2020 and place 10 for 2021
for(i in 1:8) {
  r.list2020_W[[i]] <- fc(8, i,ANPP_FORECAST_2020_W)
  r.list2020_T[[i]] <- fc(8, i,ANPP_FORECAST_2020_T)
  r.list2020_S[[i]] <- fc(8, i,ANPP_FORECAST_2020_S)
}

for(i in 1:10) {
  r.list2021_W[[i]] <- fc(10, i,ANPP_FORECAST_2021_W)
  r.list2021_T[[i]] <- fc(10, i,ANPP_FORECAST_2021_T)
  r.list2021_S[[i]] <- fc(10, i,ANPP_FORECAST_2021_S)
}

#Changes from list to DF
All_cor_2020_W <- bind_rows(r.list2020_W)
All_cor_2020_T <- bind_rows(r.list2020_T)
All_cor_2020_S <- bind_rows(r.list2020_S)

All__W <- bind_rows(r.list2021_W)
All__T <- bind_rows(r.list2021_T)
All__S <- bind_rows(r.list2021_S)

#add column with the row names
All_cor_2020_W$Forecast <- row.names(All_cor_2020_W)
All_cor_2020_T$Forecast <- row.names(All_cor_2020_T)
All_cor_2020_S$Forecast <- row.names(All_cor_2020_S)

All__W$Forecast <- row.names(All__W)
All__T$Forecast <- row.names(All__T)
All__S$Forecast <- row.names(All__S)

#Pivot to longer table and change column names

#2020
All_cor_2020_Ws <- All_cor_2020_W %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")
All_cor_2020_Ts <- All_cor_2020_T %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")
All_cor_2020_Ss <- All_cor_2020_S %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

#2021
All__Ws <- All__W %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")
All__Ts <- All__T %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")
All__Ss <- All__S %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

#Change collum names
colnames(All_cor_2020_Ws) <- c('Date','Level','Rsqure')
colnames(All_cor_2020_Ts) <- c('Date','Level','Rsqure')
colnames(All_cor_2020_Ss) <- c('Date','Level','Rsqure')

colnames(All__Ws) <- c('Date','Level','Rsqure')
colnames(All__Ts) <- c('Date','Level','Rsqure')
colnames(All__Ss) <- c('Date','Level','Rsqure')

#asDate
All_cor_2020_Ws$Date <- as.Date(All_cor_2020_Ws$Date, format="%Y-%m-%d")
All_cor_2020_Ts$Date <- as.Date(All_cor_2020_Ts$Date, format="%Y-%m-%d")
All_cor_2020_Ss$Date <- as.Date(All_cor_2020_Ss$Date, format="%Y-%m-%d")

All__Ws$Date <- as.Date(All__Ws$Date, format="%Y-%m-%d")
All__Ts$Date <- as.Date(All__Ts$Date, format="%Y-%m-%d")
All__Ss$Date <- as.Date(All__Ss$Date, format="%Y-%m-%d")

#Creates the base graphs
ggpl <- function(x){
  
    a <- ggplot(data = x, aes(x=Date)) +
    
      geom_line(aes(y=Rsqure, group=Level, colour =Level), size=0.3) +
    
      geom_point(aes(y=Rsqure, group=Level, colour =Level), size=2) +
      
      ylim(0,1) +
      
      theme_bw() +
    
      scale_x_date(date_breaks = "2 week", date_labels ="%b/%d") +
    
      scale_color_manual(values = c("33339","#CC6600","firebrick")) +
      
      labs (x = "Date of Forecast", y= "R squer of linear model") +
      
      theme(
      axis.text.x = element_text(angle = 0, vjust =0, hjust=0),
      axis.title.x =element_blank(),
      axis.title.y =element_text(face = "italic",size = 8),)
    
    return(a)
}

 #Annotated graps

# ggpl(All_cor_2020_Ws) +  labs (
#             title = "2020 Winter")
# 
# 
# ggpl(All_cor_2020_Ts) +  labs (
#             title = "2020 Transision") 
# 
# ggpl(All_cor_2020_Ss) +  labs (
#             title = "2020 SUmmer") 
# 
# 
# ggpl(All__Ws) +  labs (
#             title = "2021 Winter") 
# 
# ggpl(All__Ts) +  labs (
#             title = "2021 Transision") 
# 
# ggpl(All__Ss) +  labs (
#             title = "2021 SUmmer") 

##Summer and Spring
##Creates an empty list for the for loop
list2020_W_Spring <- list()
list2020_T_Spring <- list()
list2020_S_Spring <- list()

list2020_W_Summer <- list()
list2020_T_Summer <- list()
list2020_S_Summer <- list()

list2021_W_Spring <- list()
list2021_T_Spring <- list()
list2021_S_Spring <- list()

list2021_W_Summer <- list()
list2021_T_Summer <- list()
list2021_S_Summer <- list()

#For loop for 2020  Spring
#the final observation was on place 8 for 2020 and place 10 for 2021
for(i in 1:4) {
list2020_W_Spring[[i]] <- fc(4, i,ANPP_FORECAST_2020_W)
list2020_T_Spring[[i]] <- fc(4, i,ANPP_FORECAST_2020_T)
list2020_S_Spring[[i]] <- fc(4, i,ANPP_FORECAST_2020_S)
}

#Loop for 2020 Summer
for(i in 5:8) {
list2020_W_Summer[[i]] <- fc(8, i,ANPP_FORECAST_2020_W)
list2020_T_Summer[[i]] <- fc(8, i,ANPP_FORECAST_2020_T)
list2020_S_Summer[[i]] <- fc(8, i,ANPP_FORECAST_2020_S)
}

#For loop for 2021  Spring
for(i in 1:5) {
list2021_W_Spring[[i]] <- fc(5, i,ANPP_FORECAST_2021_W)
list2021_T_Spring[[i]] <- fc(5, i,ANPP_FORECAST_2021_T)
list2021_S_Spring[[i]] <- fc(5, i,ANPP_FORECAST_2021_S)
}

#Loop for 2021 Summer
for(i in 6:10) {
list2021_W_Summer[[i]] <- fc(10, i,ANPP_FORECAST_2021_W)
list2021_T_Summer[[i]] <- fc(10, i,ANPP_FORECAST_2021_T)
list2021_S_Summer[[i]] <- fc(10, i,ANPP_FORECAST_2021_S)
}

#Changes from list to DF 2020
All_cor_2020_W_Spring <- bind_rows(list2020_W_Spring) 
All_cor_2020_T_Spring <- bind_rows(list2020_T_Spring)
All_cor_2020_S_Spring <- bind_rows(list2020_S_Spring)

All_cor_2020_W_Summer <- bind_rows(list2020_W_Summer)
All_cor_2020_T_Summer <- bind_rows(list2020_T_Summer)
All_cor_2020_S_Summer <- bind_rows(list2020_S_Summer)

#Changes from list to DF 2021
All__W_Spring <- bind_rows(list2021_W_Spring)
All__T_Spring <- bind_rows(list2021_T_Spring)
All__S_Spring <- bind_rows(list2021_S_Spring)

All__W_Summer <- bind_rows(list2021_W_Summer)
All__T_Summer <- bind_rows(list2021_T_Summer)
All__S_Summer <- bind_rows(list2021_S_Summer)

#add column with the Forecast date
#2020
All_cor_2020_W_Spring$Forecast  <- row.names(All_cor_2020_W_Spring) 
All_cor_2020_T_Spring$Forecast  <- row.names(All_cor_2020_T_Spring)
All_cor_2020_S_Spring$Forecast  <- row.names(All_cor_2020_S_Spring)

All_cor_2020_W_Summer$Forecast  <- row.names(All_cor_2020_W_Summer)
All_cor_2020_T_Summer$Forecast  <- row.names(All_cor_2020_T_Summer)
All_cor_2020_S_Summer$Forecast  <- row.names(All_cor_2020_S_Summer)

#2021
All__W_Spring$Forecast  <- row.names(All__W_Spring)
All__T_Spring$Forecast  <- row.names(All__T_Spring)
All__S_Spring$Forecast  <- row.names(All__S_Spring)

All__W_Summer$Forecast  <- row.names(All__W_Summer)
All__T_Summer$Forecast  <- row.names(All__T_Summer)
All__S_Summer$Forecast  <- row.names(All__S_Summer)

#pivot longer and change column names
#2020 Spring
All_cor_2020_W_Spring <- All_cor_2020_W_Spring %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

All_cor_2020_T_Spring <- All_cor_2020_T_Spring %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

All_cor_2020_S_Spring <- All_cor_2020_S_Spring %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

#2020 Summer
All_cor_2020_W_Summer <- All_cor_2020_W_Summer %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

All_cor_2020_T_Summer <- All_cor_2020_T_Summer %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

All_cor_2020_S_Summer <- All_cor_2020_S_Summer %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

#2021 Spring
All__W_Spring <- All__W_Spring %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")  

All__T_Spring <- All__T_Spring %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

All__S_Spring <- All__S_Spring %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

#2021 Summer
All__W_Summer <- All__W_Summer %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

All__T_Summer <- All__T_Summer %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre") 

All__S_Summer <- All__S_Summer %>% pivot_longer(cols = c(Abv,Avg,Blw), names_to = "Time",values_to = "R_squre")

#Change the name of the colum 
colnames(All_cor_2020_W_Spring) <- c('Date','Level','Rsqure')
colnames(All_cor_2020_T_Spring) <- c('Date','Level','Rsqure')
colnames(All_cor_2020_S_Spring) <- c('Date','Level','Rsqure')

colnames(All_cor_2020_W_Summer)<- c('Date','Level','Rsqure')
colnames(All_cor_2020_T_Summer) <- c('Date','Level','Rsqure')
colnames(All_cor_2020_S_Summer) <- c('Date','Level','Rsqure')

#2021
colnames(All__W_Spring) <- c('Date','Level','Rsqure')
colnames(All__T_Spring) <- c('Date','Level','Rsqure')
colnames(All__S_Spring) <- c('Date','Level','Rsqure')

colnames(All__W_Summer)<- c('Date','Level','Rsqure')
colnames(All__T_Summer) <- c('Date','Level','Rsqure')
colnames(All__S_Summer) <- c('Date','Level','Rsqure')

#asDate
#2020
All_cor_2020_W_Spring$Date <- as.Date(All_cor_2020_W_Spring$Date, format="%Y-%m-%d")
All_cor_2020_T_Spring$Date <- as.Date(All_cor_2020_T_Spring$Date, format="%Y-%m-%d")
All_cor_2020_S_Spring$Date <- as.Date(All_cor_2020_S_Spring$Date, format="%Y-%m-%d")

All_cor_2020_W_Summer$Date <- as.Date(All_cor_2020_W_Summer$Date, format="%Y-%m-%d")
All_cor_2020_S_Spring$Date <- as.Date(All_cor_2020_T_Summer$Date, format="%Y-%m-%d")
All_cor_2020_S_Summer$Date <- as.Date(All_cor_2020_S_Summer$Date, format="%Y-%m-%d")

#2021
All__W_Spring$Date <- as.Date(All__W_Spring$Date, format="%Y-%m-%d")
All__T_Spring$Date <- as.Date(All__T_Spring$Date, format="%Y-%m-%d")
All__S_Spring$Date <- as.Date(All__S_Spring$Date, format="%Y-%m-%d")

All__W_Summer$Date <- as.Date(All__W_Summer$Date, format="%Y-%m-%d")
All__T_Summer$Date <- as.Date(All__T_Summer$Date, format="%Y-%m-%d")
All__S_Summer$Date <- as.Date(All__S_Summer$Date, format="%Y-%m-%d")


#Creates the base graphs

ggpl <- function(x){
  
    a <- ggplot(data = x, aes(x=Date)) +
    
      geom_line(aes(y=Rsqure, group=Level, colour =Level), size=0.3) +
    
      geom_point(aes(y=Rsqure, group=Level, colour =Level), size=2) +
      
      ylim(0,1) +
      
      theme_bw() +
    
      scale_x_date(date_breaks = "2 week", date_labels ="%b/%d") +
    
      scale_color_manual(values = c("33339","#CC6600","firebrick")) +
      
      labs (x = "Date of Forecast", y= "R squer of linear model") +
      
       theme_minimal() +
      
            theme(axis.title = element_blank()) 
    
    return(a)
}

#Annotated graps
SU20W <- ggpl(All_cor_2020_W_Spring) +  labs (
            title = "Winter Rain: Spring 2020") 

SP20W <- ggpl(All_cor_2020_W_Summer) +  labs (
            title = "Summer 2020") 

SU20T <- ggpl(All_cor_2020_T_Spring) +  labs (
            title = "Transition Rain: Spring 2020") 

SP20T <- ggpl(All_cor_2020_T_Summer) +  labs (
            title = "Summer 2020") 

SU20S <- ggpl(All_cor_2020_S_Spring) +  labs (
            title = "Summer Rain: Spring 2020") 

SP20S <- ggpl(All_cor_2020_S_Summer) +  labs (
            title = "Summer 2020") 


Plot0 <- SU20W + SP20W + SU20T + SP20T + SU20S + SP20S + plot_layout(ncol = 2, guides = 'collect')


#2021
#Annotated graps
SU21W <- ggpl(All__W_Spring) +  labs (
            title = "Winter Rain: Spring 2021") 

SP21W <- ggpl(All__W_Summer) +  labs (
            title = "Summer 2021") 

SU21T <- ggpl(All__T_Spring) +  labs (
            title = "Transition Rain: Spring 2021") 

SP21T <- ggpl(All__T_Summer) +  labs (
            title = "Summer 2021") 

SU21S <- ggpl(All__S_Spring) +  labs (
            title = "Summer Rain: Spring 2021") 

SP21S <- ggpl(All__S_Summer) +  labs (
            title = "Summer 2021") 


Plot1 <- SU21W + SP21W + SU21T + SP21T + SU21S + SP21S + plot_layout(ncol = 2, guides = 'collect')

Forecast Maps

## Call features from state and county
states <- map_data("state")
county <- map_data("county")

#Creates Data frames and filters states 
state_df <- states %>%  filter(str_detect(region, c("arizona","new mexico")))
county_df <- county %>%  filter(str_detect(region, c("arizona","new mexico")))

#changes state maps to simple forms to make it eay
states_sf <- sf::st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))
county_sf <- sf::st_as_sf(maps::map("county", plot = FALSE, fill = TRUE)) 

#select only AZ and NM Partial matching in the row ID
south_west_s1 <- states_sf %>% filter(str_detect(ID, c("arizona")))
south_west_s2 <- states_sf %>% filter(str_detect(ID, c("new mexico")))
south_west_s <- rbind(south_west_s1, south_west_s2)

south_west_nm <- county_sf %>% filter(str_detect(ID,"new mexico"))
south_west_az <- county_sf %>% filter(str_detect(ID,"arizona"))
south_west_c <- rbind(south_west_nm, south_west_az)

#For exprot to shhp
#south_west_c %>%
 # as_tibble() %>%
  #st_as_sf() %>% # this is the important part!
  #st_write("41421_2.shp")
######

###Create Centroids of counties to use for labels

coordinates(county_df) <- c("long", "lat")

# Get centroids
ctrs <- lapply(unique(county_df$group), 
               function(x) gCentroid(SpatialPoints(county_df[county_df$group==x,])))
ctrsout <- setNames( ctrs , unique(county_df$group ) )

##Create data frame 
centroid_df <- do.call(rbind, lapply(ctrsout, data.frame, stringsAsFactors=FALSE))

#finds and  adds unique values
uniq_logical <- !duplicated(county_df$group)
centroid_df$county <- county_df$subregion[uniq_logical]

#change names of collums
names(centroid_df) <- c("lon", "lat", 'county')

centroid_sf <-  st_as_sf(centroid_df, coords = c("lon", "lat")) %>% 
  st_set_crs(st_crs(south_west_s))

#########
#Prep data for maps

#We use 2020 year to get the lat and long
ANPP_GEOM <- subset(GC, year == '2020', select=c("gridID","lon","lat"))

#Filters Year and select columns for ANPP and drops NA values
ANPP_ALL <- subset(gfg_data, select=c("gridID", "Forecast","NPP_predict_below","NPP_predict_avg","NPP_predict_above")) %>% 
  drop_na()

jointdataset <- merge(ANPP_ALL, summerwinter, by = 'gridID')


#Slips the forecast file in to lists
ANPP_ALLFORECAST <- split(ANPP_ALL, ANPP_ALL$Forecast)

#infinite_replace <- function (x,y){
#  ANPP_ALLFORECAST[x]$y[is.infinite(ANPP_ALLFORECAST[x]$y)] <- mean(ANPP_ALLFORECAST[x]$y, na.rm = T)}

#replaces infinite values with
# infitinet_replace <- function(x) {
#   ANPP_ALLFORECAST[["8/24/2021"]][["NPP_predict_avg"]]is.infinite(ANPP_ALLFORECAST["8/24/2021"]][["NPP_predict_avg"]]) <- mean(ANPP_ALLFORECAST[["8/34/2021"]][["NPP_predict_avg"]], na.rm = T)
# }
# 
# 
# max(ANPP_ALLFORECAST[["8/24/2021"]][["NPP_predict_avg"]])



library("writexl")
write_xlsx(jointdataset,"data/ANPP_ALLFORECAST.xlsx")


##Choose Date for Forecast to draw 
#in every line drops NA values, transform df to sf object and changes the coordinates system to mach others

three_maps <- function(x){
  
#Adds coordinates to forecasts
ANPP <- merge(ANPP_GEOM, ANPP_ALLFORECAST[[x]], by = "gridID", all = TRUE) %>% drop_na() %>% st_as_sf(coords = c("lon", "lat")) %>% st_set_crs(st_crs(south_west_s))

#Draw 3 scenerios 

  y <- ggplot() +
  geom_sf() +
  geom_sf(data = ANPP, mapping = aes(col=NPP_predict_above)) +
   geom_sf(data = south_west_s, fill = NA, color=alpha("#003366",0.4), size=2.0) + 
  geom_sf(data = south_west_c, fill= NA, color=alpha("#3399FF",0.4), size=1.0) +
   #geom_label_repel(data=centroid_df, aes(lon, lat,label = county,group = county), size=2,nudge_x = -0.2, nudge_y = -0.2) +
    scale_colour_viridis_c(aesthetics = "col", option = "viridis", direction = -1, limits=c(0,2300)) +
  theme_void() +
  labs(title=paste("ANPP forecast for",str_split(names(ANPP_ALLFORECAST)[x],"/", simplify = TRUE)[3]),
      subtitle = paste(names(ANPP_ALLFORECAST)[x], "above average"), 
       col="ANNP \n (pounds/acre) \n")
  
   y <- y + theme(
      plot.title=element_text(family= "Courier",face="bold", size=18, colour = "#CC6600"),
      plot.subtitle = element_text(face="italic", size=12),
      )
  
  
  w <- ggplot() +
      geom_sf() +
      geom_sf(data = ANPP, mapping = aes(col=NPP_predict_avg)) +   scale_colour_viridis_c(aesthetics = "col", option = "viridis", direction = -1, limits=c(0,2300)) +
   geom_sf(data = south_west_s, fill = NA, color=alpha("#003366",0.4), size=2.0) + 
  geom_sf(data = south_west_c, fill= NA, color=alpha("#3399FF",0.4), size=1.0) +
  #geom_label_repel(data=centroid_df, aes(lon, lat,label = county,group = county), size=2,nudge_x = -0.2, nudge_y = -0.2) +
    
  theme_void() +
  labs(title=paste("ANPP forecast for",str_split(names(ANPP_ALLFORECAST)[x],"/", simplify = TRUE)[3]),
      subtitle = paste(names(ANPP_ALLFORECAST)[x], "average"), 
       col="ANNP \n (pounds/acre) \n")
  
  w <- w + theme(
      plot.title=element_text(family= "Courier",face="bold", size=18, colour = "#CC6600"),
      plot.subtitle = element_text(face="italic", size=12),
      )
  
  n <- ggplot() +
  geom_sf() +
  geom_sf(data = ANPP, mapping = aes(col=NPP_predict_below)) +
   geom_sf(data = south_west_s, fill = NA, color=alpha("#003366",0.4), size=2.0) + 
  geom_sf(data = south_west_c, fill= NA, color=alpha("#3399FF",0.4), size=1.0) +
  #geom_label_repel(data=centroid_df, aes(lon, lat,label = county,group = county),size=2,nudge_x = -0.2, nudge_y = -0.2) +
    scale_colour_viridis_c(aesthetics = "col", option = "viridis", direction = -1, limits=c(0,2300)) +
  theme_void() +
  labs(title=paste("ANPP forecast for",str_split(names(ANPP_ALLFORECAST)[x],"/", simplify = TRUE)[3]),
      subtitle = paste(names(ANPP_ALLFORECAST)[x], "bellow average"), 
       col="ANNP \n (pounds/acre) \n")
  
   n <- n + theme(
      plot.title=element_text(family= "Courier",face="bold", size=18, colour = "#CC6600"),
      plot.subtitle = element_text(face="italic", size=12),
      )
  
    print(y)
    print(w)
    print(n)
  
}

one_map <- function(x){
  
#Adds coordinates to forecasts
ANPP <- merge(ANPP_GEOM, ANPP_ALLFORECAST[[x]], by = "gridID", all = TRUE) %>% drop_na() %>% st_as_sf(coords = c("lon", "lat")) %>% st_set_crs(st_crs(south_west_s))

#Draw 1 scenerios 
  y <- ggplot() +
  geom_sf() +
  geom_sf(data = ANPP, mapping = aes(col=NPP_predict_above)) +
  geom_sf(data = south_west_s, fill = NA, color=alpha("#003366",0.4), size=2.0) + 
  geom_sf(data = south_west_c, fill= NA, color=alpha("#3399FF",0.4), size=1.0) +
  #geom_label_repel(data=centroid_df, aes(lon, lat,label = county,group = county),size=2,nudge_x = -0.2, nudge_y = -0.2) +
  scale_colour_viridis_c(aesthetics = "col", option = "viridis", direction = -1, limits=c(0,2300)) +
  theme_void() +
  labs(title=paste("ANPP forecast for",str_split(names(ANPP_ALLFORECAST)[x],"/", simplify = TRUE)[3]),
      subtitle = paste(names(ANPP_ALLFORECAST)[x], "final observation"), 
       col="ANNP \n (pounds/acre) \n")
  
    
   y <- y + theme(
      plot.title=element_text(family= "Courier",face="bold", size=18, colour = "#0066CC"),
      plot.subtitle = element_text(face="italic", size=12),
      )
  

  print(y)
  
}

Point Graphs

#load the files
SWC_4_14_21 <- read.csv(file = "data/SWC_4_14_21.csv", head = TRUE, sep=",")
SWC_5_15_20 <- read.csv(file = "data/SWC_5_15_20.csv", head = TRUE, sep=",")
SWC_8_24_21 <- read.csv(file = "data/SWC_8_24_21.csv", head = TRUE, sep=",")
SWC_9_12_20 <- read.csv(file = "data/SWC_9_12_20.csv", head = TRUE, sep=",")

countyaverage <- function(x) {

  date <- unique(x$Forecast)
  
  df <- summarise(group_by(x,state, county, Forecast),
                  NPP_predict_below=mean(NPP_predict_below), 
                  NPP_predict_avg=mean(NPP_predict_avg), 
                  NPP_predict_above=mean(NPP_predict_above))   
 
  df1 <- df %>% ggplot(aes(y = reorder_within(county, NPP_predict_avg, state))) +
                  geom_point(mapping = aes(x = NPP_predict_avg, 
                                           colour = 'Average')) +
                  geom_point(mapping = aes(x = NPP_predict_below,
                                           colour = 'Below')) +
                  geom_point(mapping = aes(x = NPP_predict_above,
                                           colour = 'Above')) +
                  facet_wrap(~state, scales = "free_y") +
                  scale_y_reordered() +
                  scale_color_manual(values = c("33339","#CC6600","firebrick")) +
                  labs(y = "", x="ANNP(lbs/acre)",
                  title = paste("Average NPP forecast per county",date), 
                  col="Forecast Scenarios")
print(df1)

}

Results

A dry 2020

If we can see the maps for 2020, the expectations where high, blue everywhere! But, he correlations for the first forecast where bellow 0.6 and finally the last observation turns out to be quite yellow.

This can be explains by the winter precipitations that where quite good, so the model attributed above average production of biomass for the 05/15, but this was far removed from the final forecast because of the low precipitation during the summer. Once the precipitation forecast where included in beginning of June, the model showed high correlation, that is because the scenario was stable during all the forecast season.

The point graphic shows the best and worst biomas production from each county, in the 3 different scenarios, the final observations has only green points becouse all the scenarios are e7ual

## [1] 2251.483

Moonson 2021

In 2021, the first forecast where gloomy in the three scenarios, but the last observation turn out to be very blue. The correlation coefficient for 2021 where above 0.6 at all times, but they fluctuated more than in 2021.

The forecast for 2021 showed low biomass (ANNP) at the begging of the season, because the winter rains that feed the spring growth where low.

The fluctuation in the R2 can be attributed to the monsoon of 2021. It was record breaking monsoon season and the model is not well adjusted for extremes scenarios. By the end of the summer, the rains where above normal and the correlation took a long time to adjust to the observation (because the more rain keep coming!). The above scenario was the higher correlation across all forecast in 2021.

Conclusions

There is a lot to understand about Forecast and Productivity for dry lands, as in every model, the extremes are hard to predict, the 2020 and 2021 seasons are excellent study cases to understand how far are the forecast from unexpected situations. The rains for winter are easier to predict than the monsoon season (Smith et al., 2019), and the model gets more accurate as the monsoon season forecast improves or shows below normal precipitation from the start.

References

Brown, J. R., & Thorpe, J. (2008). Climate Change and Rangelands: Responding Rationally to Uncertainty. 4.

Hartman, M. D., Parton, W. J., Derner, J. D., Schulte, D. K., Smith, W. K., Peck, D. E., Day, K. A., Del Grosso, S. J., Lutz, S., Fuchs, B. A., Chen, M., & Gao, W. (2020). Seasonal grassland productivity forecast for the U.S. Great Plains using Grass‐Cast. Ecosphere, 11(11). https://doi.org/10.1002/ecs2.3280

ILRI. (2021). Rangelands ATLAS. International Livestock Research Institute. https://www.rangelandsdata.org/atlas/sites/default/files/2021-05/Rangelands%20Atlas.pdf

National Drought Mitigation Center. 2022. [Data Set] available online: https://grasscast.unl.edu/Archive.aspx

NOAA. (2021). GHCN (Global Historical Climatology Network) [Data Set] available online: https://www.ncdc.noaa.gov/cdo-web/

Seager, R., Ting, M., Held, I., Kushnir, Y., Lu, J., Vecchi, G., Huang, H.-P., Harnik, N., Leetmaa, A., Lau, N.-C., Li, C., Velez, J., & Naik, N. (2007). Model Projections of an Imminent Transition to a More Arid Climate in Southwestern North America. Science, 316(5828), 1181–1184. https://doi.org/10.1126/science.1139601

Smith, W. K., Dannenberg, M. P., & Yang, D. (2019). Remote sensing of dryland ecosystem structure and function_ Progress, challenges, and opportunities. Remote Sensing of Environment, 23.