In danamelamed/swnsmodelr: Model weather station data in Southwest Nova Scotia

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)

Appendix B: Graphs

Model Stations

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('')

    }
}

Validation Stations

#
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('')

    }
}