library(swnsmodelr) library(ggplot2) knitr::opts_chunk$set(eval = FALSE)
# get evi values # Make raster df tables evi_raster_df <- make_temporal_raster_df(file.path("C:",'MScAG',"EVI_200m"), start_date = ymd('2012-04-01'), end_date = ymd('2012-11-30'), date_chars = c(5,-4), date_format = "%Y_%j") %>% add_date_columns() %>% filter(month >= 4 & month <= 11) gdd_raster_df <- make_temporal_raster_df(file.path("F:","output",'gdd10','monthly'), start_date = ymd('2012-04-01'), end_date = ymd('2012-11-30'), date_chars = c(7,-5), date_format = "%Y_%m_%d")
For each monthly gdd raster... plot the corresponding EVI raster
for(i in seq_along(gdd_raster_df[[1]])){ gdd_raster_row <- filter(gdd_raster_df,date_time == gdd_raster_df[[2]][[i]]) gdd_raster <- raster(gdd_raster_row[[1]]) evi_raster_row <- filter(evi_raster_df,date_time == evi_raster_df[[2]][[i]]) evi_raster <- raster(evi_raster_row[[1]]) plot(evi_raster) plot(gdd_raster) }
evi_df <- extract_temporal_raster_values(evi_raster_df, df_12, col_name = "evi", verbose = TRUE) %>% dplyr::select(stationid, EASTING, NORTHING, date_time, evi) %>% filter(date_time %in% evi_raster_df$date_time) %>% mutate(evi_calib = (evi + 2000)/12000*1000)
# add gdd values gdd_evi_df <- extract_temporal_raster_values(gdd_raster_df, evi_df, col_name = "gdd", verbose = TRUE)
Add GDD zone values
gdd_zones <- raster(file.path("F:","data","rasters","swns_ecodist.tif")) names(gdd_zones) <- 'gdd_zones' gdd_evi_df <- extract_constant_raster_values(gdd_evi_df, gdd_zones)
gdd_zones_list <- gdd_zones %>% unique() for(j in seq_along(gdd_zones_list)){ gdd_zone_now <- gdd_zones_list[[j]] gdd_evi_df_now <- gdd_evi_df %>% filter(gdd_zones == gdd_zone_now) stations_list <- gdd_evi_df_now %>% filter(!duplicated(stationid)) %>% dplyr::select(stationid) stations_list <- stations_list$stationid gdd_zone_now_raster <- gdd_zones gdd_zone_now_raster[gdd_zone_now_raster != gdd_zone_now] <- 0 for(i in seq_along(stations_list)){ print('') gdd_evi_df_now <- gdd_evi_df %>% filter(gdd_zones == gdd_zone_now) gdd_evi_df_now <- gdd_evi_df_now %>% filter(stationid == stations_list[[i]]) if(length(gdd_evi_df_now[[1]] != 0)){ # plot point on map point_now <- gdd_evi_df_now coordinates(point_now) = ~ EASTING + NORTHING plot(gdd_zone_now_raster, main = paste0("Station : ",stations_list[[i]]," (Zone ",gdd_zone_now,")"), col = c("grey","pink"), legend = FALSE) points(point_now, pch = 21, bg = "black", lwd = 1.5) print(ggplot(data = gdd_evi_df_now) + geom_point(aes(x = yday, y= gdd, colour = "blue")) + geom_smooth(aes(x = yday, y=gdd, colour = "blue"), se = FALSE) + geom_point(aes(x = yday, y= evi_calib, colour = "green")) + geom_smooth(aes(x = yday, y=evi_calib, colour = "green"), se = FALSE) + facet_wrap(~ year, ncol = 2) + scale_colour_discrete(labels = c("GDD0","EVI")) + guides(colour = guide_legend(title = "Variable")) + labs(x = "Day of the Year", y = "Value") + theme(text = element_text(size=12)) ) }else{} print('') print('') } }
# gdd_zones_list <- list(1,2,3,4,5) for(j in seq_along(gdd_zones_list)){ gdd_zone_now <- gdd_zones_list[[j]] gdd_evi_df_now <- gdd_evi_df %>% filter(gdd_zones == gdd_zone_now & stationid %in% val_stations_list) stations_list <- gdd_evi_df_now %>% filter(!duplicated(stationid)) %>% dplyr::select(stationid) stations_list <- stations_list$stationid print(stations_list) gdd_zone_now_raster <- gdd_zones gdd_zone_now_raster[gdd_zone_now_raster != gdd_zone_now] <- 0 for(i in seq_along(stations_list)){ print('') gdd_evi_df_now <- gdd_evi_df %>% filter(gdd_zones == gdd_zone_now & stationid %in% val_stations_list) gdd_evi_df_now <- gdd_evi_df_now %>% filter(stationid == stations_list[[i]]) if(length(gdd_evi_df_now[[1]] != 0)){ # plot point on map point_now <- gdd_evi_df_now coordinates(point_now) = ~ EASTING + NORTHING plot(gdd_zone_now_raster, main = paste0("Station : ",stations_list[[i]]," (Zone ",gdd_zone_now,")"), col = c("grey","pink"), legend = FALSE) points(point_now, pch = 21, bg = "black", lwd = 1.5) print(ggplot(data = gdd_evi_df_now) + geom_point(aes(x = yday, y= gdd, colour = "blue")) + geom_smooth(aes(x = yday, y=gdd, colour = "blue"), se = FALSE) + geom_point(aes(x = yday, y= evi_calib, colour = "green")) + geom_smooth(aes(x = yday, y=evi_calib, colour = "green"), se = FALSE) + facet_wrap(~ year, ncol = 2) + scale_colour_discrete(labels = c("GDD0","EVI")) + guides(colour = guide_legend(title = "Variable")) + labs(x = "Day of the Year", y = "Value") + theme(text = element_text(size=12)) ) }else{} print('') print('') } }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.