## ----------------------------------------------------------------------------- # 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)))) ## ----------------------------------------------------------------------------- # 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),] ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- # 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" ## ----------------------------------------------------------------------------- # 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 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" ## ----------------------------------------------------------------------------- # 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" ## ----------------------------------------------------------------------------- # 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')