## ----setup_1, include = FALSE------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width=6, fig.height=4 ) options(scipen = 9999) ## ----libs, message=FALSE, warning=FALSE--------------------------------------- library(sf) library(dplyr) library(ncdfgeom) ## ----data_shape_2------------------------------------------------------------- prcp_data <- readRDS(system.file("extdata/climdiv-pcpndv.rds", package = "ncdfgeom")) print(prcp_data, max_extra_cols = 0) plot(prcp_data$date, prcp_data$`0101`, col = "red", xlab = "date", ylab = "monthly precip (inches)", main = "Sample Timeseries for 0101-'Northern Valley'") lines(prcp_data$date, prcp_data$`0101`) ## ----data_shape--------------------------------------------------------------- climdiv_poly <- read_sf(system.file("extdata/climdiv.gpkg", package = "ncdfgeom")) print(climdiv_poly) plot(st_geometry(climdiv_poly), main = "Climate Divisions with 0101-'Northern Valley' Highlighted") plot(st_geometry(filter(climdiv_poly, CLIMDIV == "0101")), col = "red", add = TRUE) ## ----write_ts, warning = FALSE------------------------------------------------ climdiv_centroids <- climdiv_poly %>% st_transform(5070) %>% # Albers Equal Area st_set_agr("constant") %>% st_centroid() %>% st_transform(4269) %>% #NAD83 Lat/Lon st_coordinates() %>% as.data.frame() nc_file <- "climdiv_prcp.nc" prcp_dates <- prcp_data$date prcp_data <- select(prcp_data, -date) prcp_meta <- list(name = "climdiv_prcp_inches", long_name = "Estimated Monthly Precipitation (Inches)") write_timeseries_dsg(nc_file = nc_file, instance_names = climdiv_poly$CLIMDIV, lats = climdiv_centroids$Y, lons = climdiv_centroids$X, times = prcp_dates, data = prcp_data, data_unit = rep("inches", (ncol(prcp_data) - 1)), data_prec = "float", data_metadata = prcp_meta, attributes = list(title = "Demonstation of ncdfgeom"), add_to_existing = FALSE) climdiv_poly <- st_sf(st_cast(climdiv_poly, "MULTIPOLYGON")) write_geometry(nc_file = "climdiv_prcp.nc", geom_data = climdiv_poly, variables = "climdiv_prcp_inches") ## ----ncdump------------------------------------------------------------------- try({ncdump <- system(paste("ncdump -h", nc_file), intern = TRUE) cat(ncdump, sep = "\n")}, silent = TRUE) ## ----read--------------------------------------------------------------------- # First read the timeseries. prcp_data <- read_timeseries_dsg("climdiv_prcp.nc") # Now read the geometry. climdiv_poly <- read_geometry("climdiv_prcp.nc") ## ----data_summary------------------------------------------------------------- names(prcp_data) class(prcp_data$time) names(prcp_data$varmeta$climdiv_prcp_inches) prcp_data$data_unit prcp_data$data_prec str(names(prcp_data$data_frames$climdiv_prcp_inches)) prcp_data$global_attributes names(climdiv_poly) ## ----p_colors_source, echo=FALSE---------------------------------------------- # Because we've gotta have pretty colors! p_colors <- function (n, name = c("precip_colors")) { # Thanks! https://quantdev.ssri.psu.edu/tutorials/generating-custom-color-palette-function-r p_rgb <- col2rgb(c("#FAFBF3", "#F0F8E3", "#D4E9CA", "#BBE0CE", "#B7DAD0", "#B0CCD7", "#A9B8D7", "#A297C2", "#8F6F9E", "#684A77", "#41234D")) precip_colors = rgb(p_rgb[1,],p_rgb[2,],p_rgb[3,],maxColorValue = 255) name = match.arg(name) orig = eval(parse(text = name)) rgb = t(col2rgb(orig)) temp = matrix(NA, ncol = 3, nrow = n) x = seq(0, 1, , length(orig)) xg = seq(0, 1, , n) for (k in 1:3) { hold = spline(x, rgb[, k], n = n)$y hold[hold < 0] = 0 hold[hold > 255] = 255 temp[, k] = round(hold) } palette = rgb(temp[, 1], temp[, 2], temp[, 3], maxColorValue = 255) palette } ## ----plot, fig.height=6, fig.width=8------------------------------------------ climdiv_poly <- climdiv_poly %>% st_transform(3857) %>% # web mercator st_simplify(dTolerance = 5000) title <- paste0("\n Sum of: ", prcp_data$varmeta$climdiv_prcp_inches$long_name, "\n", format(prcp_data$time[1], "%Y-%m", tz = "UTC"), " - ", format(prcp_data$time[length(prcp_data$time)], "%Y-%m", tz = "UTC")) prcp_sum <- apply(prcp_data$data_frames$climdiv_prcp_inches, 2, sum, na.rm = TRUE) prcp <- data.frame(CLIMDIV = names(prcp_sum), prcp = as.numeric(prcp_sum), stringsAsFactors = FALSE) %>% right_join(climdiv_poly, by = "CLIMDIV") %>% st_as_sf() plot(prcp["prcp"], lwd = 0.1, pal = p_colors, breaks = seq(0, 14000, 1000), main = title, key.pos = 3, key.length = lcm(20)) ## ----setup_dontrun, eval = FALSE---------------------------------------------- # # Description here: ftp://ftp.ncdc.noaa.gov/pub/data/cirs/climdiv/divisional-readme.txt # prcp_url <- "ftp://ftp.ncdc.noaa.gov/pub/data/cirs/climdiv/climdiv-pcpndv-v1.0.0-20190408" # # prcp_file <- "prcp.txt" # # download.file(url = prcp_url, destfile = prcp_file, quiet = TRUE) # # division_url <- "ftp://ftp.ncdc.noaa.gov/pub/data/cirs/climdiv/CONUS_CLIMATE_DIVISIONS.shp.zip" # division_file <- "CONUS_CLIMATE_DIVISIONS.shp.zip" # # download.file(url = division_url, destfile = division_file, quiet = TRUE) # zip::unzip("CONUS_CLIMATE_DIVISIONS.shp.zip") # # climdiv_poly <- read_sf("GIS.OFFICIAL_CLIM_DIVISIONS.shp") %>% # select("CLIMDIV", CLIMDIV_NAME = "NAME") %>% # mutate(CLIMDIV = ifelse(nchar(as.character(CLIMDIV)) == 3, # paste0("0",as.character(CLIMDIV)), # as.character(CLIMDIV))) %>% # st_simplify(dTolerance = 0.0125) # # month_strings <- c("jan", "feb", "mar", "apr", "may", "jun", # "jul", "aug", "sep", "oct", "nov", "dec") # # prcp_data <- read.table(prcp_file, header = FALSE, # colClasses = c("character", # rep("numeric", 12)), # col.names = c("meta", month_strings)) # # # Here we gather the data into a long format and prep it for ncdfgeom. # prcp_data <- prcp_data %>% # gather(key = "month", value = "precip_inches", -meta) %>% # mutate(climdiv = paste0(substr(meta, 1, 2), substr(meta, 3, 4)), # year = substr(meta, 7, 10), # precip_inches = ifelse(precip_inches < 0, NA, precip_inches)) %>% # mutate(date = as.Date(paste0("01", "-", month, "-", year), # format = "%d-%b-%Y")) %>% # select(-meta, -month, -year) %>% # filter(climdiv %in% climdiv_poly$CLIMDIV) %>% # spread(key = "climdiv", value = "precip_inches") %>% # filter(!is.na(`0101`)) # as_tibble() # # # Now make sure things are in the same order. # climdiv_names <- names(prcp_data)[2:length(names(prcp_data))] # climdiv_row_order <- match(climdiv_names, climdiv_poly$CLIMDIV) # climdiv_poly <- climdiv_poly[climdiv_row_order, ] # # sf::write_sf(climdiv_poly, "climdiv.gpkg") # saveRDS(prcp_data, "climdiv-pcpndv.rds") # # unlink("GIS*") # unlink("CONUS_CLIMATE_DIVISIONS.shp.zip") # unlink("prcp.txt") ## ----p_colors, eval=FALSE----------------------------------------------------- # p_colors <- function (n, name = c("precip_colors")) { # # Thanks! https://quantdev.ssri.psu.edu/tutorials/generating-custom-color-palette-function-r # p_rgb <- col2rgb(c("#FAFBF3", "#F0F8E3", "#D4E9CA", # "#BBE0CE", "#B7DAD0", "#B0CCD7", # "#A9B8D7", "#A297C2", "#8F6F9E", # "#684A77", "#41234D")) # precip_colors = rgb(p_rgb[1,],p_rgb[2,],p_rgb[3,],maxColorValue = 255) # name = match.arg(name) # orig = eval(parse(text = name)) # rgb = t(col2rgb(orig)) # temp = matrix(NA, ncol = 3, nrow = n) # x = seq(0, 1, , length(orig)) # xg = seq(0, 1, , n) # for (k in 1:3) { # hold = spline(x, rgb[, k], n = n)$y # hold[hold < 0] = 0 # hold[hold > 255] = 255 # temp[, k] = round(hold) # } # palette = rgb(temp[, 1], temp[, 2], temp[, 3], maxColorValue = 255) # palette # } ## ----cleanup, echo=FALSE------------------------------------------------------ unlink("climdiv_prcp.nc")