knitr::opts_chunk$set(eval = FALSE) library(swnsmodelr) library(knitr)
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)
# 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 } }
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.