diff --git a/.github/workflows/secretScan.yml b/.github/workflows/secretScan.yml new file mode 100644 index 0000000..1af40d6 --- /dev/null +++ b/.github/workflows/secretScan.yml @@ -0,0 +1,13 @@ +name: gitleaks + +on: [push, pull_request] + +jobs: + gitleaks: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + with: + fetch-depth: '0' + - name: gitleaks-action + uses: gitleaks/gitleaks-action@v1.6.0 diff --git a/DESCRIPTION b/DESCRIPTION index 0bac300..04a7f8b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,9 +1,10 @@ Package: coldpool Type: Package Title: AFSC/RACE Groundfish Assessment Program EBS and NBS temperature products -Version: 2.1 +Version: 2.2-1 Authors@R: c(person("Sean", "Rohan", email = "sean.rohan@noaa.gov", role = c("aut", "cre")), person("Lewis", "Barnett", email = "lewis.barnett@noaa.gov", role = c("aut", "ctb")), + person("Kelly", "Kearney", role = c("ctb")), person("Emily", "Markowitz", role = c("ctb"))) Maintainer: Sean Rohan Depends: R (>= 3.5), @@ -28,7 +29,7 @@ Description: This package calculates the area of the cold pool in the eastern Be License: GPL (>=2) Encoding: UTF-8 LazyData: false -RoxygenNote: 7.2.1 +RoxygenNote: 7.2.3 Imports: colorspace, viridis Suggests: knitr, rmarkdown, RODBC, getPass, readr, testthat, cowplot, metR, cowplot diff --git a/NAMESPACE b/NAMESPACE index edb1759..c1424dc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(calcindices_temp) export(cold_pool_index) export(compare_cpa_station_filter) export(convert_ddm_to_dd) @@ -21,6 +22,7 @@ export(make_raster_stack) export(nbs_ebs_bottom_temperature) export(nbs_ebs_surface_temperature) export(nbs_mean_temperature) +export(parseregion) export(plot_loocv_rmse) export(scale_fill_fermenter_custom) export(sql_to_rqry) diff --git a/R/data2raster.R b/R/data2raster.R new file mode 100644 index 0000000..0c65cf3 --- /dev/null +++ b/R/data2raster.R @@ -0,0 +1,437 @@ +#' Interpolate data points to a masked raster +#' +#' This function is a variant of the interpolate_variable function that allows +#' for user-set masking polygons and choice of fitting method. +#' +#' @param dat Input data frame. Must include columns for latitude, longitude, +#' and variable to be interpolated. +#' @param dat.year Year +#' @param var.col Character vector denoting the name of the variable column. +#' @param lat.col Character vector indicating the name of the latitude column. +#' @param lon.col Character vector indicating the name of the longitude column. +#' @param in.crs Character vector (PROJ4 or WKT2/PROJ6 string) indicating the +#' coordinate reference system for the data. +#' @param interpolation.crs Character vector (PROJ4 or WKT2/PROJ6 string) +#' indicating the coordinate reference system to use for the interpolation. +#' Historically EPSG:3338 +#' @param cell.resolution Dimensions of interpolation grid cell in meters. +#' @param nm Maximum number of cells to use for interpolation. +#' @param select.region Region for interpolation as a character string. Can be +#' either a string accepted by akgfmaps::get_base_layers or a simple feature +#' polygon +#' @param fittypes vector of strings specifying which fit algorithms to use: +#' "nn" = nearest neighbor +#' "idwnmax4" = inverse distance weighting with nmax=4 +#' "idw" = inverse distance weighting +#' "exp" = ordinary kriging w/ exponential VGM, +#' "sph" = ordinary kriging w/ spherical VGM +#' "bes" = ordinary kriging w/ bessel VGM +#' "gau" = ordinary kriging w/ gaussian VGM +#' "cir" = ordinary kriging w/ circular VGM +#' "mat" = ordinary kriging w/ matern VGM +#' "ste" = ordinary kriging w/ Stein's matern VGM +#' "tps" = thin plate spline +#' @param bbox Bounding box for the interpolated grid. Can be either a string +#' ("sebs", "bs.south", "nebs", "ebs", "bs.north", "bs.all") corresponding to +#' specific bounds, or a 4-element named vector with values xmin, xmax, ymin, +#' and ymax. +#' @return raster stack of interpolated data, with one layer per input fit type + +data2raster <- function(dat, + dat.year, + var.col, + lat.col, + lon.col, + in.crs = "+proj=longlat +datum=NAD83", + interpolation.crs, + cell.resolution, + nm = Inf, + select.region = "sebs", + fittypes = c("nn", "idwnmax4", "idw", "exp", "sph", "bes", "gau", "cir", "mat", "ste", "tps"), + bbox = NULL) +{ + + # Rename data columns for ease + + names(dat)[which(names(dat) == var.col)] <- "var.col" + names(dat)[which(names(dat) == lat.col)] <- "lat.col" + names(dat)[which(names(dat) == lon.col)] <- "lon.col" + + # Remove NAs + + if(any(is.na(dat$var.col))) { + print(paste0("Removing ", + sum(is.na(dat$var.col)), + " var.col NA values from data set")) + dat <- dat %>% + dplyr::filter(!is.na(var.col)) + } + + # Load AKGF polygon for masking, or the user-supplied mask + + mymask <- parseregion(select.region, interpolation.crs) + + # Set dimensions for raster cells --- + # Note: If SEBS or EBS/NEBS is specified, we set the raster to match the older + # ArcGIS rasters. If user does not specify anything, a raster is created that + # encompasses the masking polygon plus a bit of a buffer and aligns with the + # default SEBS raster. + + if (is.character(bbox)) { + if (bbox %in% c("sebs","bs.south")) { + bbox = c(xmin = -1625000, + xmax = -35000, + ymin = 379500, + ymax = 1969500) + } else if (bbox %in% c("nebs","ebs","bs.north","bs.all")) { + extrap.box <- c(xmn = -179.5, xmx = -157, ymn = 50, ymx = 68) + plot.boundary <- akgfmaps::transform_data_frame_crs(data.frame(x = c(extrap.box['xmn'], extrap.box['xmx']), + y = c(extrap.box['ymn'], extrap.box['ymx'])), + out.crs = interpolation.crs) + bbox = c(xmin = plot.boundary$x[1], + xmax = plot.boundary$x[2], + ymin = plot.boundary$y[1], + ymax = plot.boundary$y[2], + ) + } else { + stop("Unrecognized value of bbox") + } + } else if (is.null(bbox)) { + bbox <- st_bbox(st_buffer(mymask, cell.resolution*5)) + bbox <- round(bbox/cell_resolution)*cell_resolution # round to align rows/cols + xshift <- ((-1625000/cell_resolution) - floor(-1625000/cell_resolution))*cell_resolution + yshift <- ((379500/cell_resolution) - floor(379500/cell_resolution))*cell_resolution + bbox = bbox + c(xshift,yshift,xshift,yshift) + } + + sp_interp.raster <- raster::raster(xmn = bbox["xmin"], + xmx = bbox["xmax"], + ymn = bbox["ymin"], + ymx = bbox["ymax"], + nrows = floor((bbox["ymax"]-bbox["ymin"])/cell.resolution), + ncols = floor((bbox["xmax"]-bbox["xmin"])/cell.resolution)) + + raster::projection(sp_interp.raster) <- interpolation.crs + + # Transform data for interpolation ---- + sp_interp.df <- unique(dat) + sp::coordinates(sp_interp.df) <- c(x = "lon.col", y = "lat.col") + sp::proj4string(sp_interp.df) <- sp::CRS(in.crs) + sp_interp.df <- sp::spTransform(sp_interp.df, sp::CRS(interpolation.crs)) + + + # Fit data and build rasters + + # IDW base for variogram-based fits + + if (any(c("exp", "sph", "bes", "gau", "cir", "mat", "ste") %in% fittypes)){ + idw_vgm_fit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + nmax = Inf) + } + + # Loop over requested fit types + + outraster <- stack() + + for (f in fittypes) { + + # Fit data + + if (f == "nn") { + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + set = list(idp = 0), + nmax = 4) + } + if (f == "idwnmax4") { + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + set = list(idp = 2), + nmax = 4) + } + if (f == "idw") { + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + set = list(idp = 2), + nmax = Inf) + } + if (f == "exp") { + vgfit <- gstat::fit.variogram(variogram(idw_vgm_fit), + vgm(c("Exp"))) + + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + model = vgfit, + nmax = nm) + } + if (f == "sph") { + vgfit <- gstat::fit.variogram(variogram(idw_vgm_fit), + vgm(c("Sph"))) + + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + model = vgfit, + nmax = nm) + } + if (f == "bes") { + vgfit <- gstat::fit.variogram(variogram(idw_vgm_fit), + vgm(c("Bes"))) + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + model = vgfit, + nmax = nm) + } + if (f == "gau") { + vgfit <- gstat::fit.variogram(variogram(idw_vgm_fit), + vgm(c("Gau"))) + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + model = vgfit, + nmax = nm) + } + if (f == "cir") { + vgfit <- gstat::fit.variogram(variogram(idw_vgm_fit), + vgm(c("Cir"))) + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + model = vgfit, + nmax = nm) + } + if (f == "mat") { + + vgfit <- gstat::fit.variogram(variogram(idw_vgm_fit), + vgm(c("Mat"))) + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + model = vgfit, + nmax = nm) + } + if (f == "ste") { + + vgfit <- gstat::fit.variogram(variogram(idw_vgm_fit), + vgm(c("Ste"))) + + xfit <- gstat::gstat(formula = var.col ~ 1, + locations = sp_interp.df, + model = vgfit, + nmax = nm) + } + if (f == "tps") { + xfit <- fields::Tps(sp::coordinates(sp_interp.df), + sp_interp.df$var.col) + } + + # Build predictor + + if (f == "tps") { + xpredict <- raster::interpolate(sp_interp.raster, + xfit) + } else { + xpredict <- predict(xfit, as(sp_interp.raster, "SpatialGrid")) + } + + # Interpolate data to raster and mask outside of polygon + + rastmp <- xpredict %>% + akgfmaps::rasterize_and_mask(mymask) + + outraster <- addLayer(outraster, rastmp) + } + return(outraster) +} + + + +#' Parse region string or simple feature +#' +#' This helper function parses whether region is provided as a string or as a +#' simple feature. If the former, it loads the appropriate simple feature using +#' the akgfmaps::get_base_layers function; if the latter, it re-projects the +#' provided features to the specified coordinate reference system. +# +#' @param nameorpoly string specifying base layer name, or simple feature object +#' @param crs coordinate reference system to apply to simple feature object +#' @return simple feature object in the requested CRS +#' @export + +parseregion <- function(nameorpoly, crs) { + + if (is.character(nameorpoly)) { + region_layers <- akgfmaps::get_base_layers(select.region = nameorpoly, + set.crs = crs) + mymask <- region_layers$survey.area + + } else if (is.data.frame(nameorpoly)) { + + mymask <- sf::st_transform(nameorpoly, crs = sf::st_crs(crs)) + } else { + stop("region must be an akgfmaps-compatible string or shapefile-derived dataframe") + } + + return(mymask) + +} + + + +#' Calculate temperature-related indices for the eastern Bering Sea shelf +#' +#' This routine replicates the primary mean temperature and area-below-x indices +#' as found in the cold_pool_index data frame, but allows for user-specified +#' variations. +#' +#' @param datafile Full filename of .csv file holding point-sampled surface and +#' bottom temperature data +#' @param maskpoly String specifying masking polygon from +#' akgfmaps::get_base_layers, or simple feature object to use as masking +#' polygon +#' @param select_years vector of years for which indices will be calculated +#' @param cell_resolution resolution (m) to use for interpolation raster +#' @param method fitting method to use for interpolation (see data2raster +#' fittype parameter for options) +#' @param proj_crs CRS to use for interpolation, or NULL to use coldpool:::ebs_proj_crs +#' @param threshold vector of temperature threshold values (deg C) for which +#' area-less-than-x indices will be calculated +#' @param bbox Bounding box to use for interpolated raster (see data2raster) +#' @param bottomvar Name of bottom temperaure variable column in the datafile, +#' or NULL to skip calculation of bottom temperature metrics +#' @param surfacevar Name of surface temperaure variable column in the datafile, +#' or NULL to skip calculation of surface temperature metrics +#' @return data frame with the following column +#' YEAR: year of indices +#' MEAN_GEAR_TEMPERATURE: mean interpolated and masked bottomvar +#' MEAN_SURFACE_TEMPERATURE: mean interpolated and masked surfacevar +#' MEAN_BT_LT100M: mean interpolated bottomvar within the intersection of the +#' <100m survey strata and masked region +#' AREA_LTE_X: total area (km^2) where temperature is less than X deg C +#' @export + +calcindices_temp <- function (datafile, + maskpoly, + select_years = NULL, + cell_resolution = 5000, + method = "ste", + proj_crs = NULL, + threshold=c(-1,0,1,2), + bbox = NULL, + bottomvar = "gear_temperature", + surfacevar = "surface_temperature") +{ + # Read temperature data from .csv file + + print("Reading temperature data...") + + temperature_df <- read.csv(file = datafile, + stringsAsFactors = FALSE) + + # Replace NULL input with defaults + + year_vec <- sort(unique(temperature_df$year)) + if(!is.null(select_years)) { + year_vec <- year_vec[year_vec %in% select_years] + } + + + if (is.null(proj_crs)) { + proj_crs <- coldpool:::ebs_proj_crs + } + + # Mask polygons: main and <100m polygon + + print("Checking/building masking polygon") + + maskpoly <- parseregion(maskpoly, proj_crs) + + ebs_layers <- akgfmaps::get_base_layers(select.region = "sebs", set.crs = proj_crs) + + lt100_strata <- ebs_layers$survey.strata %>% + dplyr::filter(Stratum %in% c(10, 20, 31, 32, 41, 42, 43)) %>% + dplyr::group_by(SURVEY) %>% + dplyr::summarise() + + lt100_strata <- sf::st_intersection(lt100_strata, maskpoly) + + # Initialize data frame + + print("Initializing data frame") + + nyr <- length(select_years) + nthresh <- length(threshold) + + area_lte <- matrix(NA, nyr, nthresh) + + bt_df <- data.frame(YEAR = numeric(length = nyr), + MEAN_GEAR_TEMPERATURE = numeric(length = nyr), + MEAN_BT_LT100M = numeric(length = nyr), + MEAN_SURFACE_TEMPERATURE = numeric(length = nyr)) + + # Build rasters for surface and gear temperature + + + for (i in 1:nyr) { + + print(paste("Calculating indices:", select_years[i])) + + bt_df$YEAR[i] <- select_years[i] + + # Bottom temperature + + if (!is.null(bottomvar)) { + + ras <- data2raster(dat = dplyr::filter(temperature_df, year == select_years[i]), + dat.year = yr, + in.crs = "+proj=longlat", + interpolation.crs = proj_crs, + cell.resolution = cell_resolution, + lon.col = "longitude", + lat.col = "latitude", + var.col = bottomvar, + nm = Inf, + select.region = maskpoly, + fittypes = method, + bbox = bbox) + ras <- ras[[1]] # 1-layer rasterstack to plain raster + + bt_df$MEAN_GEAR_TEMPERATURE[i] <- raster::values(ras) %>% + mean(na.rm = TRUE) + + lt100_temp <- raster::mask(ras, lt100_strata) + bt_df$MEAN_BT_LT100M[i] <- mean(lt100_temp@data@values, na.rm = TRUE) + + for (it in 1:nthresh) { + area_lte[i,it] <- ras %>% + cpa_from_raster(raster_units = "m", temperature_threshold = threshold[it]) + } + + } + + # Surface temperature + + if (!is.null(surfacevar)) { + ras <- data2raster(dat = dplyr::filter(temperature_df, year == select_years[i]), + dat.year = yr, + in.crs = "+proj=longlat", + interpolation.crs = proj_crs, + cell.resolution = cell_resolution, + lon.col = "longitude", + lat.col = "latitude", + var.col = surfacevar, + nm = Inf, + select.region = maskpoly, + fittypes = method, + bbox = bbox) + ras <- ras[[1]] # 1-layer rasterstack to plain raster + + bt_df$MEAN_SURFACE_TEMPERATURE[i] <- raster::values(ras) %>% + mean(na.rm = TRUE) + } + } + area_lte <- as.data.frame(area_lte) + names(area_lte) <- gsub("-", "m", sprintf("AREA_LTE_%.1f", threshold)) + + bt_df <- cbind(bt_df, area_lte) + + return(bt_df) + +} \ No newline at end of file diff --git a/R/sysdata.rda b/R/sysdata.rda index 94c3d88..c5ced6e 100644 Binary files a/R/sysdata.rda and b/R/sysdata.rda differ diff --git a/analysis/roms_level3_coldpool_part2.R b/analysis/roms_level3_coldpool_part2.R new file mode 100644 index 0000000..774cda8 --- /dev/null +++ b/analysis/roms_level3_coldpool_part2.R @@ -0,0 +1,107 @@ +# roms_level3_coldpool_part2.R +# +# This script is part of the post-processing suite for the Bering10K ROMS +# hindcast simulation. It calculates the survey-based and survey-replicated +# cold pool indices and adds those to the existing cold pool index file +# (created by roms_level3_coldpool_part1.m) + +library(coldpool) +library(ncdf4) +library(lubridate) +# source('extended_utils.R') Functions are now included in the package + +# Index file + +moxdir = "~/Documents/Research/Working/mox_bumblereem/" +fname <- file.path(moxdir, "roms_for_public", "B10K-K20_CORECFS", "Level3", "B10K-K20_CORECFS_coldpool.nc") +# fname = here::here("test.nc") +ncin <- nc_open(fname) + +# Read existing data and parse years we need to calculate + +time <- ncvar_get(ncin,"time") +tdate <- lubridate::ymd("1900-01-01") + lubridate::days(time) + +meanbtmp <- ncvar_get(ncin,"average_bottom_temp") +cpindex <- ncvar_get(ncin, "cold_pool_index") +thresh <- ncvar_get(ncin, "threshold") + +nc_close(ncin) + +ismiss <- meanbtmp[1,2,] > 1e36 # Don't know how to query fill value in R... + +select_years <- year(tdate[ismiss]) +select_years <- select_years[select_years >= 1982 & select_years != 2020] +# select_years <- c(2005, 2017, 2018, 2019) + +# User-set variables + +proj_crs <- coldpool:::ebs_proj_crs +cell_resolution <- 5000 + +# Masks: w/ and w/o NW (82+90), w/ and w/o narrowed-to-model-200m-isobath + +akSEBS <- akgfmaps::get_base_layers(select.region = "bs.south", set.crs = proj_crs) + +b10_lte_200 <- sf::st_read(here::here("data", "B10K_lte200m_polygon.shp"), quiet = TRUE) %>% + st_set_crs("+proj=longlat") %>% + sf::st_transform(crs = sf::st_crs(proj_crs)) + +akSEBSnonw <- akSEBS$survey.strata %>% + dplyr::filter(Stratum <= 62) %>% + dplyr::group_by(SURVEY) %>% + dplyr::summarise() + +mymask <- list(akSEBSnonw, + akSEBS$survey.area, + st_intersection(akSEBSnonw, b10_lte_200), + st_intersection(akSEBS$survey.area, b10_lte_200)) + +# Calculate indices using survey and survey-replicated data + +datafile <- here::here("data", "survey_replicates_B10K-K20_CORECFS.csv") + +cpname <- gsub("-", "m", sprintf("AREA_LTE_%.1f", thresh)) + +yrmask <- year(tdate) %in% select_years +for (ii in 1:length(mymask)) { + + totarea <- as.numeric(st_area(mymask[[ii]]))/1e6 # m^2 to km^2 + + svy <- calcindices_temp(datafile, mymask[[ii]], + select_years = select_years, + bottomvar="gear_temperature", + surfacevar=NULL, + threshold=thresh) + + srep <- calcindices_temp(datafile, mymask[[ii]], + select_years = select_years, + bottomvar="model_bottom_temp", + surfacevar=NULL, + threshold=thresh) + + meanbtmp[ii,2,yrmask] <- svy$MEAN_GEAR_TEMPERATURE + meanbtmp[ii,3,yrmask] <- srep$MEAN_GEAR_TEMPERATURE + + cpindex[,ii,2,yrmask] <- aperm(data.matrix(svy[,cpname]), c(2,1))/totarea + cpindex[,ii,3,yrmask] <- aperm(data.matrix(srep[,cpname]), c(2,1))/totarea +} + +# Write to file + +ncout <- nc_open(fname, write=TRUE) + +ncvar_put(ncout, "average_bottom_temp", meanbtmp) +ncvar_put(ncout, "cold_pool_index", cpindex) + +hisstr <- ncatt_get(ncout, 0, "history") +if (length(select_years) == 1) { + newhis <- paste0(date(), ": Survey and survey-replicated data added for year ", select_years) +} else { + newhis <- paste0(date(), ": Survey and survey-replicated data added for year ", min(select_years), "-", max(select_years)) +} +hisstr <- paste(newhis, hisstr$value, sep="\n") + +ncatt_put(ncout, 0, "history", hisstr) + +nc_close(ncout) diff --git a/man/calcindices_temp.Rd b/man/calcindices_temp.Rd new file mode 100644 index 0000000..ba3f200 --- /dev/null +++ b/man/calcindices_temp.Rd @@ -0,0 +1,61 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data2raster.R +\name{calcindices_temp} +\alias{calcindices_temp} +\title{Calculate temperature-related indices for the eastern Bering Sea shelf} +\usage{ +calcindices_temp( + datafile, + maskpoly, + select_years = NULL, + cell_resolution = 5000, + method = "ste", + proj_crs = NULL, + threshold = c(-1, 0, 1, 2), + bbox = NULL, + bottomvar = "gear_temperature", + surfacevar = "surface_temperature" +) +} +\arguments{ +\item{datafile}{Full filename of .csv file holding point-sampled surface and +bottom temperature data} + +\item{maskpoly}{String specifying masking polygon from +akgfmaps::get_base_layers, or simple feature object to use as masking +polygon} + +\item{select_years}{vector of years for which indices will be calculated} + +\item{cell_resolution}{resolution (m) to use for interpolation raster} + +\item{method}{fitting method to use for interpolation (see data2raster +fittype parameter for options)} + +\item{proj_crs}{CRS to use for interpolation, or NULL to use coldpool:::ebs_proj_crs} + +\item{threshold}{vector of temperature threshold values (deg C) for which +area-less-than-x indices will be calculated} + +\item{bbox}{Bounding box to use for interpolated raster (see data2raster)} + +\item{bottomvar}{Name of bottom temperaure variable column in the datafile, +or NULL to skip calculation of bottom temperature metrics} + +\item{surfacevar}{Name of surface temperaure variable column in the datafile, +or NULL to skip calculation of surface temperature metrics} +} +\value{ +data frame with the following column + YEAR: year of indices + MEAN_GEAR_TEMPERATURE: mean interpolated and masked bottomvar + MEAN_SURFACE_TEMPERATURE: mean interpolated and masked surfacevar + MEAN_BT_LT100M: mean interpolated bottomvar within the intersection of the + <100m survey strata and masked region + AREA_LTE_X: total area (km^2) where temperature is less than X deg C +} +\description{ +This routine replicates the primary mean temperature and area-below-x indices +as found in the cold_pool_index data frame, but allows for user-specified +variations. +} diff --git a/man/data2raster.Rd b/man/data2raster.Rd new file mode 100644 index 0000000..e9e08f5 --- /dev/null +++ b/man/data2raster.Rd @@ -0,0 +1,74 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data2raster.R +\name{data2raster} +\alias{data2raster} +\title{Interpolate data points to a masked raster} +\usage{ +data2raster( + dat, + dat.year, + var.col, + lat.col, + lon.col, + in.crs = "+proj=longlat +datum=NAD83", + interpolation.crs, + cell.resolution, + nm = Inf, + select.region = "sebs", + fittypes = c("nn", "idwnmax4", "idw", "exp", "sph", "bes", "gau", "cir", "mat", "ste", + "tps"), + bbox = NULL +) +} +\arguments{ +\item{dat}{Input data frame. Must include columns for latitude, longitude, +and variable to be interpolated.} + +\item{dat.year}{Year} + +\item{var.col}{Character vector denoting the name of the variable column.} + +\item{lat.col}{Character vector indicating the name of the latitude column.} + +\item{lon.col}{Character vector indicating the name of the longitude column.} + +\item{in.crs}{Character vector (PROJ4 or WKT2/PROJ6 string) indicating the +coordinate reference system for the data.} + +\item{interpolation.crs}{Character vector (PROJ4 or WKT2/PROJ6 string) +indicating the coordinate reference system to use for the interpolation. +Historically EPSG:3338} + +\item{cell.resolution}{Dimensions of interpolation grid cell in meters.} + +\item{nm}{Maximum number of cells to use for interpolation.} + +\item{select.region}{Region for interpolation as a character string. Can be +either a string accepted by akgfmaps::get_base_layers or a simple feature +polygon} + +\item{fittypes}{vector of strings specifying which fit algorithms to use: +"nn" = nearest neighbor +"idwnmax4" = inverse distance weighting with nmax=4 +"idw" = inverse distance weighting +"exp" = ordinary kriging w/ exponential VGM, +"sph" = ordinary kriging w/ spherical VGM +"bes" = ordinary kriging w/ bessel VGM +"gau" = ordinary kriging w/ gaussian VGM +"cir" = ordinary kriging w/ circular VGM +"mat" = ordinary kriging w/ matern VGM +"ste" = ordinary kriging w/ Stein's matern VGM +"tps" = thin plate spline} + +\item{bbox}{Bounding box for the interpolated grid. Can be either a string +("sebs", "bs.south", "nebs", "ebs", "bs.north", "bs.all") corresponding to +specific bounds, or a 4-element named vector with values xmin, xmax, ymin, +and ymax.} +} +\value{ +raster stack of interpolated data, with one layer per input fit type +} +\description{ +This function is a variant of the interpolate_variable function that allows +for user-set masking polygons and choice of fitting method. +} diff --git a/man/parseregion.Rd b/man/parseregion.Rd new file mode 100644 index 0000000..aa9f92d --- /dev/null +++ b/man/parseregion.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data2raster.R +\name{parseregion} +\alias{parseregion} +\title{Parse region string or simple feature} +\usage{ +parseregion(nameorpoly, crs) +} +\arguments{ +\item{nameorpoly}{string specifying base layer name, or simple feature object} + +\item{crs}{coordinate reference system to apply to simple feature object} +} +\value{ +simple feature object in the requested CRS +} +\description{ +This helper function parses whether region is provided as a string or as a +simple feature. If the former, it loads the appropriate simple feature using +the akgfmaps::get_base_layers function; if the latter, it re-projects the +provided features to the specified coordinate reference system. +} diff --git a/plots/2022_stacked_temperature_simple_label.png b/plots/2022_stacked_temperature_simple_label.png index 21ca167..55df4da 100644 Binary files a/plots/2022_stacked_temperature_simple_label.png and b/plots/2022_stacked_temperature_simple_label.png differ