Retrieve data from GYGA¶

João Vasco Silva, CIMMYT-Zimbabwe
Marloes van Loon, WUR

Introduction¶

This notebook complements an earlier notebook describing the methodology for yield gap decomposition. That earlier notebook makes use of water-limited yield data to decompose yield gaps. Such data were derived using the scripts documented in this notebook. The reader is referred to that earlier notebook for further information about the concepts and definitions considered in yield gap analysis. To make the approach fully reproducible, it is explained here how to retrieve the water-limited yield data from the Global Yield Gap Atlas (GYGA) using available APIs for acquiring such data. Also here an example is provided for for wheat in Ethiopia.

Load required R packages¶

First, the R packages needed to run this workflow are loaded.

# package names
packages <- c("dplyr", "tidyr", "httr", "jsonlite", "sf", "reshape2")
#
# install packages
installed_packages <- packages %in% rownames(installed.packages())
if(any(installed_packages == FALSE)){
  install.packages(packages[!installed_packages], repos="http://cran.us.r-project.org", quiet=T)}
#
# load packages
invisible(lapply(packages, function(x) suppressMessages(require(x, character.only=T, quietly=T, warn.conflicts=F))))

Access to GYGA data¶

An example on how data can be accessed from GYGA, and linked to farmer field data, using an API specifically set for wheat water-limited yields in Ethiopia. The chunk of code below illustrates how to access from the GYGA API the following data for a given country: (a) crop yield data for all weather stations, (b) weather station coordinates and information and, (c) climate zones. Please note an internet connection is needed to access the API.

# retrieve data for all weather stations
request <- httr::GET(
  "https://www.yieldgap.org/apigyga/json/cropcountrystationyear?accesstoken=anonymous&par1=30&par2=5",
  httr::add_headers())
wheat_data <- httr::content(request, as = "text")
wheat_data <- jsonlite::fromJSON(wheat_data)[["items"]]
wheat_data$harvestyear <- paste0('Yw_', wheat_data$harvestyear)
wheat_data <- dcast(wheat_data, country + country_id + climatezone +
                                station + station_id + crop + crop_id ~ harvestyear,
                    value.var='yw')
#
# get coordinates of weather stations
request <- httr::GET("https://www.yieldgap.org/apigyga/metadata/stations",
                     httr::add_headers())
json_data <- httr::content(request, as="text")
gyga_stations <- jsonlite::fromJSON(json_data)[["items"]]
gyga_stations <- gyga_stations[gyga_stations$station_id %in% unique(wheat_data$station_id),]
#
# retrieve climate zones gyga
request <- httr::GET("https://www.yieldgap.org/apigyga/metadata/climatezones",
                     httr::add_headers())
json_data <- httr::content(request, as="text")
cz_data <- jsonlite::fromJSON(json_data)[["items"]]
cz_data <- cz_data[cz_data$country == unique(wheat_data$country),]

The next step entails merging the yield data from GYGA with the respective weather station coordinates and climate zone information. Once data are merged, the crop yield data is made spatially explicit using the st_as_sf() function. Please refer to the chunk of code below.

# merge wheat data with gyga stations
wheat_data <- merge(wheat_data[,c("country", "country_id", "station", "station_id",
                                  "climatezone", "crop", "crop_id",
                                  "Yw_1998", "Yw_1999", "Yw_2000", "Yw_2001",
                                  "Yw_2002", "Yw_2003", "Yw_2004", "Yw_2005",
                                  "Yw_2006", "Yw_2007", "Yw_2008", "Yw_2009",
                                  "Yw_2010", "Yw_2011", "Yw_2012", "Yw_2013",
                                  "Yw_2014", "Yw_2015", "Yw_2016", "Yw_2017")],
                    gyga_stations[,c("station_id", "latitude", "longitude")],
                    by="station_id", all.x=TRUE)
#
# merge wheat data with climate zones
wheat_data <- merge(wheat_data,
                    cz_data[,c("climatezone", "climatezone_id")],
                    by="climatezone", all.x=TRUE)
#
# convert to spatial dataframe
wheat_data <- st_as_sf(wheat_data, coords=c("longitude", "latitude")) %>%
  st_set_crs(4326)

Link to farmer field data¶

Once the GYGA data are retrieved, it is then possible to merge for each field available in the farmer field data. To do so, the data frame with the farmer field data is loaded and made spatially explicit based on the latitude and longitude values of each field. The climate zones for the country of interest are also loaded and added to the farmer field data using the st_join() function. This is described in the chunk of code below.

# read .csv file with data
file <- 'https://raw.githubusercontent.com/jvasco323/EiA_YGD_workflow/main/data-gps-coordinates-original.csv'
data <- read.csv(url(file))
data <- unique(data)
#
# convert df to spatial dataframe
data <- data[!is.na(data$LON),]
data <- data[!is.na(data$LAT),]
data <- st_as_sf(data, coords=c("LON", "LAT")) %>%
  st_set_crs(4326)
#
# load climate zones and merge to field data
file <- 'https://raw.githubusercontent.com/jvasco323/EiA_YGD_workflow/main/data-climate-zone-eth.rds'
climate_zone <- readRDS(url(file))
#
# join field data with climate zones
data <- st_join(data, climate_zone, join=st_intersects, left=FALSE)

The next step is to get the closest weather station to each field within the climate zone of the field. The closest weather station must be within the same climate zone of the field, but when climate zones are very fragmented, it is advisable to consider the closest weather station independently of the climate zone (more information on this is provided below). The chunk of code below adds the closest weather station to each field.

# get the closest station within each climate zone
data_rfwh  <- NULL
for(i in unique(data$GYGA_CZ)){
  #
  # subset each climate zone
  dat <- subset(data, GYGA_CZ==i)
  dat_grid <- subset(wheat_data, climatezone==i)
  #
  # find the nearest station within each climate zone
  find_nearest_RWS_CZ <- st_join(dat, dat_grid, join=st_nearest_feature, left=FALSE)
  nearest <- st_nearest_feature(dat, dat_grid)
  dist <- st_distance(dat, dat_grid[nearest,], by_element=TRUE)
  dist <- data.frame(dist)
  #
  # conver m to km
  dist$dist <- as.numeric(dist$dist)/1000 #convert from m to km
  find_nearest_RWS_CZ <- cbind(find_nearest_RWS_CZ, dist)
  #
  # bind all data from the loop
  data_rfwh <- rbind(data_rfwh, find_nearest_RWS_CZ)
}
# add -99 to yw not available
data_rfwh[c(22:41)][is.na(data_rfwh[c(22:41)])] <- -99
## Warning in `[<-.data.frame`(`*tmp*`, c(22:41), value = structure(list(Yw_1998 =
## c(-99, : provided 21 variables to replace 20 variables

Now, the fields with water-limited yield data from the same climate zone of the weather station are identified. To do this, a subset is created containing the farmer field data for each the water-limited yield is greater than 0. A new column is added to that data frame to flag the source of the water-limited yield data.

# fields with yield data
data_rfwh_with_data <- subset(data_rfwh, Yw_2010 >= 0) # arbitrary year of the data
data_rfwh_with_data$data_from <- "same_cz"

Some fields do not fall inside climate zones for which crop model simulations were done and hence, they do not have a weather station coupled. These fields are identified through a subset of the data set containing the fields for which the water-limited yield is below 0 (i.e., equal to -99). For these fields, the distance from each field to the closest weather station is estimated using the chunk of code below.

# fields with no yield data
data_rfwh_no_data <- subset(data_rfwh, Yw_2010 < 0) # arbitrary year of the data
data_rfwh_no_data <- data_rfwh_no_data[,c('hhid')]
#
# check GYGA station closer than 30km
nearest <- st_nearest_feature(data_rfwh_no_data, wheat_data)
dist <- st_distance(data_rfwh_no_data, wheat_data[nearest,], by_element=TRUE)
dist <- data.frame(dist)
dist$dist <- as.numeric(dist$dist)/1000 #convert from m to km
data_rfwh_no_data <- cbind(data_rfwh_no_data, dist)

Fields with a weather station closest than 30km are identified first. For those fields, the yield data from the respective weather station are used. Also here a new column is added to the data frame to flag the source of the water-limited yield data. This is done in the chunk of code below.

# fields closer than 30km: get data from respective station
data_rfwh_no_data_30 <- subset(data_rfwh_no_data, dist <= 30)
data_rfwh_no_data_30 <- st_join(data_rfwh_no_data_30, wheat_data, join=st_nearest_feature,
                                left=FALSE)
data_rfwh_no_data_30$data_from <- "cz_station"

Lastly, crop yield data need to be retrieved for the fields located more than 30km away from any weather station. The national average is used for these fields. The chunk of code below first subsets the fields with no weather station closer than 30km, then retrieves the national average for the country of interest using the GYGA API and, finally, merges that data with the field data.

# fields further than 30km: use country average
data_rfwh_no_data_m30 <- subset(data_rfwh_no_data, dist > 30)
data_rfwh_no_data_m30$country <- 'Ethiopia'
#
# retrieve wheat data Ethiopia
request <- httr::GET(
  "https://www.yieldgap.org/apigyga/json/cropcountryyear?accesstoken=anonymous&par1=30&par2=5",
  httr::add_headers())
wheat_data_eth <- httr::content(request, as="text")
wheat_data_eth <- jsonlite::fromJSON(wheat_data_eth)[["items"]]
wheat_data_eth$harvestyear <- paste0('Yw_', wheat_data_eth$harvestyear)
wheat_data_eth <- dcast(wheat_data_eth, country + country_id +
                                        crop + crop_id ~ harvestyear, value.var='yw')
#
# join yw data of whole country with fields which have no data
data_rfwh_no_data_m30 <- merge(data_rfwh_no_data_m30,
                               wheat_data_eth[,c("country", "country_id", "crop", "crop_id",
                                                 "Yw_1998", "Yw_1999", "Yw_2000", "Yw_2001",
                                                 "Yw_2002", "Yw_2003", "Yw_2004", "Yw_2005",
                                                 "Yw_2006", "Yw_2007", "Yw_2008", "Yw_2009",
                                                 "Yw_2010", "Yw_2011", "Yw_2012", "Yw_2013",
                                                 "Yw_2014", "Yw_2015", "Yw_2016", "Yw_2017")],
                               by="country", all.x=TRUE)
data_rfwh_no_data_m30$data_from <- "country_average"

Export file with GYGA data¶

The last step in this workflow is to bring back the different subsets of fields into a single data frame. Recall: the different subsets of fields are as follows: (a) fields with a weather station within the respective climate zone, (b) fields without a climate zone for which crop model simulations were done but with a weather station closer than 30km and, (c) fields without a climate zone for which crop model simulations were done and with no weather station closer than 30km. This is implemented in the chunk of code below, where this final data frame is saved to disk as a csv file.

# final data frame
data_rfwh_final <- rbind(data_rfwh_with_data[,c("hhid", "country", "crop",
                                                "data_from", "geometry",
                                                "Yw_1998", "Yw_1999", "Yw_2000", "Yw_2001",
                                                "Yw_2002", "Yw_2003", "Yw_2004", "Yw_2005",
                                                "Yw_2006", "Yw_2007", "Yw_2008", "Yw_2009",
                                                "Yw_2010", "Yw_2011", "Yw_2012", "Yw_2013",
                                                "Yw_2014", "Yw_2015", "Yw_2016", "Yw_2017")],
                         data_rfwh_no_data_30[,c("hhid", "country", "crop",
                                                 "data_from", "geometry",
                                                 "Yw_1998", "Yw_1999", "Yw_2000", "Yw_2001",
                                                 "Yw_2002", "Yw_2003", "Yw_2004", "Yw_2005",
                                                 "Yw_2006", "Yw_2007", "Yw_2008", "Yw_2009",
                                                 "Yw_2010", "Yw_2011", "Yw_2012", "Yw_2013",
                                                 "Yw_2014", "Yw_2015", "Yw_2016", "Yw_2017")],
                         data_rfwh_no_data_m30[,c("hhid", "country", "crop",
                                                  "data_from", "geometry",
                                                  "Yw_1998", "Yw_1999", "Yw_2000", "Yw_2001",
                                                  "Yw_2002", "Yw_2003", "Yw_2004", "Yw_2005",
                                                  "Yw_2006", "Yw_2007", "Yw_2008", "Yw_2009",
                                                  "Yw_2010", "Yw_2011", "Yw_2012", "Yw_2013",
                                                  "Yw_2014", "Yw_2015", "Yw_2016", "Yw_2017")])
#
# get cimate zones for each field
data_rfwh_final <- st_join(data_rfwh_final, climate_zone, join=st_intersects, left=FALSE)
#
# remove geometry
data_rfwh_final <- st_drop_geometry(data_rfwh_final)
#
# calculate average Yw
data_rfwh_final$Yw_average <- rowMeans(data_rfwh_final[c(5:24)], na.rm=TRUE)
#
# save csv file
write.csv(data_rfwh_final, 'data-gps-coordinates-final.csv')

Final remarks¶

Please note this script does not add data to fields without GPS coordinates reported. The national average could also be considered for those fields, but this is not implemented in this script.

Data from GYGA can be made available to the OneCGIAR initiative on Excellence in Agronomy on a demand basis. Requests should be made to prof. Martin van Ittersum (martin.vanittersum@wur.nl). The script documented here is reproducible and can be used to retrieved GYGA data for other crop x country combinations. Any questions or suggestions for improving the scripts presented in this document should be addressed to j.silva@cgiar.org or marloes.vanloon@wur.nl.