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

knitr::opts_chunk$set(eval = FALSE)
library(swnsmodelr)
library(knitr)

The Input Data

Polygon Boundaries

ecodists_shp <- readOGR("F:\\data","swns_ecodist")
dem500m <- raster('f:\\data\\rasters\\input_rasters\\dem500m.tif')
#ecodist_raster <- rasterize(ecodists_shp,dem500m, field=1)
# gdd_dists <- raster(file.path("E:","Daily","gdd_dists.tif"))
# plot(gdd_dists,
#      col = rev(rainbow(n = 5,start = 0,end = 0.75)),
#      lwd = 0.1,
#      main = "Polygon boundaries: GDD5 dists")
# 
# plot(dists_sp, add = TRUE)

Rasters

# Names of rasters in Rasters directory (everything before extension)
# rasters_names_list <- list("dem",  # elevation
#                           "ptoc", # proxmity to coast
#                           "east", # easting
#                           "north",# northing
#                           "asp",  # aspect
#                            "tpi", # topographic position index
#                           "slope") # slope
# # Make into rasters objects list
# rasters_list <- lapply(FUN = raster,
#                        X   = paste0("E:\\Packages\\swnsmodelr\\Rasters\\200\\",rasters_names_list, ".tif"))
# 
# # Plot brick of rasters
# rasters_list %>% brick() %>% plot()

Make list of output rasters to analyze

raster_summary <- list() 

gdd10_rasters_names_in <- list.files('F:\\output\\500\\apr_nov',full.names = TRUE,pattern='gdd10_20')
gdd10_rasters <- list()
gdd10_rasters_names_out <- list()
c <- 1
for(i in seq_along(gdd10_rasters_names_in)){
  if(substring(gdd10_rasters_names_in[[i]],nchar(gdd10_rasters_names_in[[i]])-3,nchar(gdd10_rasters_names_in[[i]]))=='.tif'){
    print(c)
    gdd10_rasters[[c]] <- gdd10_rasters_names_in[[i]] %>% raster
    gdd10_rasters_names_out[[c]] <- gdd10_rasters_names_in[[i]] 
    c <- c +1
  }
}

# Create a list of ecodistricts
dists <- ecodists_shp@data$ECODISTRIC %>% unique()

# Create a summary data.frame to populate
dists_summary_df <- data.frame(matrix(NA, 
                                      nrow = length(dists)*length(gdd10_rasters),
                                      ncol = 5))
names(dists_summary_df) <- c('dist','mean','sd','n','raster_names')

gdd10_rasters_names_out %>% 
  unlist %>%
  substring(nchar(gdd10_rasters_names_out)-13,nchar(gdd10_rasters_names_out))


c <- 1
for(j in seq_along(gdd10_rasters)){

  raster_now <-gdd10_rasters[[j]]
  raster_name <- names(raster_now)
  print(raster_name)
   for(i in  seq_along(dists)){

      # subset dists 
      dist_now <- subset(ecodists_shp, ecodists_shp@data$ECODISTRIC == dists[[i]])
      print(dists[[i]])
      # extract gdd10 values to district
      values <- raster::extract(raster_now,
                                dist_now) %>% unlist()
      dists_summary_df$raster_names[[c]] <- raster_name
      dists_summary_df$dist[[c]] <- dists[i]
      dists_summary_df$mean[[c]] <- mean(values,na.rm=TRUE)
      dists_summary_df$sd[[c]] <- sd(values,na.rm=TRUE)
      dists_summary_df$n[[c]] <- length(values)
      c <- c+1
    }

}

Summary table

dists_summary_df <- data.frame(matrix(NA, 
                                      nrow = length(dists),
                                      ncol = length(rasters_list)))
names(dists_summary_df) <- lapply(FUN = names, X = rasters_list)
row.names(dists_summary_df) <- dists

for(i in seq_along(dists_summary_df)){ # cols (rasters)
  for(j in seq_along(dists_summary_df[[1]])){ # rows (dists)

    pixels_now <- dists_summary[[i]][[j]]$values
    mean_now <- pixels_now %>% mean(na.rm=TRUE) %>% as.numeric() %>% round(2)
    sd_now   <- pixels_now %>% sd(na.rm=TRUE) %>% as.numeric() %>% round(2)

    dists_summary_df[[i]][[j]] <- paste0(mean_now, " (",sd_now,")")

  }
}
knitr::kable(dists_summary_df,
      row.names = dists)

dists_summary_df <- data.frame(matrix(NA, 
                                      nrow = length(dists),
                                      ncol = length(rasters_list)))
names(dists_summary_df) <- names(rasters_list)
row.names(dists_summary_df) <- dists

for(i in seq_along(rasters_list)){ # cols (rasters)
  for(j in seq_along(dists)){ # rows (dists)

    pixels_now <- dists_summary[[i]][[j]]$values
    mean_now <- mean(pixels_now)
    sd_now   <- sd(pixels_now)

    dists_summary_df[j, i] <- paste0(round(mean_now, 2), " (", round(sd_now, 2), ")")

  }
}

knitr::kable(dists_summary_df,
             row.names = dists)