USGS-R/netcdf.dsg
'NetCDF' Geometry and Time Series

Package index
Search the USGS-R/netcdf.dsg package
Vignettes
Functions
64
Source code
22
Man pages
8
Home
/
GitHub
/
USGS-R/netcdf.dsg
/
R/write_geometry.R

R/write_geometry.R
In USGS-R/netcdf.dsg: 'NetCDF' Geometry and Time Series

Defines functions put_hole_node_count put_part_node_count put_node_count put_poly write_geom_POLYGON put_line write_geom_LINESTRING write_geom_POINT put_meta_nc init_ncdf get_type get_crs write_geom_data.sfc_MULTIPOLYGON write_geom_data.sfc_POLYGON write_geom_data.sfc_MULTILINESTRING write_geom_data.sfc_LINESTRING write_geom_data.sfc_POINT write_geom_data write_geometry

Documented in write_geometry 

#' @title Write geometries and attributes to NetCDF-CF
#' @param nc_file \code{character} file path to the nc file to be created.
#' @param geom_data sf \code{data.frame} with POINT, LINESTRING, MULTILINESTRING, 
#' POLYGON, or MULTIPOLYGON geometries. Note that three dimensional geometries 
#' are not supported.
#' @param instance_dim_name \code{character} Not required if adding geometry to a 
#' NetCDF-CF Discrete Sampling Geometries timeSeries file. For a new file, will 
#' use package default -- "instance" -- if not supplied.
#' @param variables \code{character} If a an existing netCDF files is provided, this 
#' list of variables that should be related to the geometries.
#'
#' @description
#' Creates a file with point, line or polygon instance data ready for the 
#' extended NetCDF-CF timeSeries featuretype format.
#' 
#' Will also add attributes if provided data has them.
#'
#'@references
#' \enumerate{
#'   \item \url{http://cfconventions.org/cf-conventions/cf-conventions.html}
#'}
#'
#' @importFrom sf st_set_geometry st_geometry
#' @export
#' 
#' @examples 
#'
#' hucPolygons <- sf::read_sf(system.file('extdata','example_huc_eta.json', package = 'ncdfgeom'))
#'
#' hucPolygons_nc <- ncdfgeom::write_geometry(nc_file=tempfile(), 
#'                                            geom_data = hucPolygons)
#' try({
#'   ncdump <- system(paste("ncdump -h", hucPolygons_nc), intern = TRUE)
#'   cat(ncdump ,sep = "\n")
#' }, silent = TRUE)
#' 
write_geometry = function(nc_file, geom_data, instance_dim_name = NULL, variables = list()) {

  type <- get_type(geom_data)
    
  if(grepl("MULTIPOINT", type)) {
  	stop("Multi point not supported yet.")
  } else if(!grepl("POLYGON", type) & 
            !grepl("LINESTRING", type) &
            type != "POINT") {
    stop("Did not find supported spatial data.")
  }
  
  att_data <- st_set_geometry(geom_data, NULL)
  geom_data <- st_geometry(geom_data)

  if(is.null(instance_dim_name)) {
    instance_dim_name <- pkg.env$instance_dim_name
  }

	if(ncol(att_data) > 0) {
		itemp <- sapply(att_data, is.factor)
		att_data[itemp] <- lapply(att_data[itemp], as.character)
		
		nc_file <- write_attribute_data(nc_file, att_data, instance_dim_name)
		
		variables <- c(variables, names(att_data))
	}

  return(write_geom_data(geom_data, nc_file, instance_dim_name, variables = variables))
}

#'@title Put geometry data in a NetCDF-CF File
#'@param nc_file A string file path to the nc file to be created. It must already have
#'an instance dimension.
#'@param geom_data An object of class \code{SpatialLines}, \code{SpatialPolygons}
#'or their sf with WGS84 lon in the x coordinate and lat in the y coordinate.
#'Note that three dimensional geometries are not supported.
#'@param instance_dim_name A string to name the instance dimension.  Defaults to "instance"
#'@param variables A character vector of variable names that the geometry data
#'container variable name will be added to.

#'@description
#'Creates a file with point, line or polygon instance data ready for the extended NetCDF-CF timeSeries featuretype format.
#'Will also add attributes if a sp dataframe object is passed in.
#'
#'@references
#'https://github.com/twhiteaker/netCDF-CF-simple-geometry
#'
#' @importFrom RNetCDF open.nc close.nc create.nc var.put.nc att.put.nc
#' @importFrom stats setNames
#' @importFrom sf st_geometry_type st_crs st_coordinates st_zm
#' @noRd
write_geom_data <- function(geom_data, nc_file, instance_dim_name, variables = c())
  UseMethod("write_geom_data")

#' @noRd
#' @name write_geom_data
write_geom_data.sfc_POINT <- function(geom_data, nc_file, instance_dim_name, variables = c()) {
  crs <- get_crs(geom_data)
  
  geom_data <- st_coordinates(geom_data)
  
  node_count <- c(1:nrow(geom_data))
  
  x_vals <- geom_data[,1]
  y_vals <- geom_data[,2]
  
  node_dim_name <- instance_dim_name
  
  nc <- init_ncdf(nc_file, x_vals, y_vals, instance_dim_name)
  
  write_geom_POINT(nc)
  
  put_meta_nc(nc, crs, variables)
  
  close.nc(nc)
  
  return(nc_file)
}

#' @noRd
#' @name write_geom_data
write_geom_data.sfc_LINESTRING <- function(geom_data, nc_file, instance_dim_name, variables = c()) {
  crs <- get_crs(geom_data)
  
  geom_data <- st_coordinates(geom_data)
  
  n_geoms <- length(unique(geom_data[, 3]))
  
  node_count <- rep(NA, n_geoms)
  x_vals <- rep(NA, nrow(geom_data))
  y_vals <- x_vals
  coord <- 1
  
  for(geom in 1:n_geoms) {
    coords <- geom_data[geom_data[, 3] == geom, ][, c(1,2)]
    nrc <- nrow(coords)

    x_vals[coord:(coord + nrc - 1)] <- coords[,1]
    y_vals[coord:(coord + nrc - 1)] <- coords[,2]
    
    node_count[geom] <- nrc
    coord <- coord + nrc
  }
  
  nc <- init_ncdf(nc_file, x_vals, y_vals, pkg.env$node_dim_name)
  
  write_geom_LINESTRING(nc, node_count, instance_dim_name, NULL)
  
  put_meta_nc(nc, crs, variables)
  
  close.nc(nc)
  
  return(nc_file)
}

#' @noRd
#' @name write_geom_data
write_geom_data.sfc_MULTILINESTRING <- function(geom_data, nc_file, 
                                                instance_dim_name, variables = c()) {
  if(grepl("Z|M", class(st_geometry(geom_data)[[1]])[1])) {
    warning("Found more than two dimensions in geometry. Removing Z and M content.")
    geom_data <- st_zm(geom_data)
  }
  
  crs <- get_crs(geom_data)
  
  geom_data <- st_coordinates(geom_data)
  
  part_dif <- diff(geom_data[, 3])
  part_dif[part_dif < 0] <- 0
  geom_dif <- diff(geom_data[, 4])
  geom_data[, 3] <- c(1, cumsum(part_dif + geom_dif) + 1)
  
  n_geoms <- length(unique(geom_data[, 4]))
  n_parts <- length(unique(geom_data[, 3]))
  
  nCoords <- nrow(geom_data)
  
  part_node_count <- rep(NA, n_parts)
  node_count <- rep(NA, n_geoms)
  
  x_vals <- rep(NA, nCoords)
  y_vals <- x_vals
  nc_part <- 0
  coord <- 1
  
  for(geom in 1:n_geoms) {
    nd_count <- 0
    
    g_data <- geom_data[geom_data[, 4] == geom, ]

    for(g_part in unique(g_data[, 3])) {
      nc_part <- nc_part + 1
      
      coords <- g_data[g_data[, 3] == g_part, c(1,2)]
      nrc <- nrow(coords)
      
      x_vals[coord:(coord+nrc-1)] <- coords[,1]
      y_vals[coord:(coord+nrc-1)] <- coords[,2]
      
      nd_count <- nd_count + nrc
      part_node_count[nc_part] <- nrc
      coord <- coord + nrc
    }
    node_count[geom] <- nd_count
  }
  
  nc <- init_ncdf(nc_file, x_vals, y_vals, pkg.env$node_dim_name)
  
  write_geom_LINESTRING(nc, node_count, instance_dim_name, part_node_count)
  
  put_meta_nc(nc, crs, variables)
  
  close.nc(nc)
  
  return(nc_file)
}

#' @noRd
#' @name write_geom_data
write_geom_data.sfc_POLYGON <- function(geom_data, nc_file, instance_dim_name, variables = c()) {
  crs <- get_crs(geom_data)
  
  geom_data <- st_coordinates(geom_data)
  
  hole_dif <- diff(geom_data[, 3])
  hole_dif[hole_dif < 0] <- 0
  
  n_geoms <- length(unique(geom_data[, 4]))
  n_parts <- n_geoms + sum(hole_dif)
  
  nCoords <- nrow(geom_data)
  
  part_node_count <- rep(NA, n_parts)  
  part_type <- rep(pkg.env$multi_val, n_parts)
  node_count <- rep(NA, n_geoms)
  
  x_vals <- rep(NA, nCoords)
  y_vals <- x_vals
  nc_part <- 0
  coord <- 1
  
  for(geom in 1:n_geoms) {
    nd_count <- 0
    g_data <- geom_data[geom_data[, 4] == geom, ]
    
    nc_parts <- length(unique(g_data[, 3]))
    
    for(g_part in 1:nc_parts) {
      nc_part <- nc_part + 1
      
      if(g_part > 1) {
        part_type[nc_part] <- pkg.env$hole_val
      }
      
      coords <- g_data[g_data[, 3] == g_part, c(1,2)]
      nrc <- nrow(coords)
      
      x_vals[coord:(coord+nrow(coords)-1)]<-coords[1:nrow(coords), 1]
      y_vals[coord:(coord+nrow(coords)-1)]<-coords[1:nrow(coords), 2]
        
      nd_count <- nd_count + nrc
      part_node_count[nc_part] <-  nrc
      coord <- coord + nrc
    }
    node_count[geom] <- nd_count
  }
  
  nc <- init_ncdf(nc_file, x_vals, y_vals, pkg.env$node_dim_name)
  
  if(all(part_type == pkg.env$multi_val)) part_node_count <- part_type <- NULL
  
  write_geom_POLYGON(nc, node_count, instance_dim_name, part_node_count, part_type)
  
  put_meta_nc(nc, crs, variables)
  
  close.nc(nc)
  
  return(nc_file)
}

#' @noRd
#' @name write_geom_data
write_geom_data.sfc_MULTIPOLYGON <- function(geom_data, nc_file, instance_dim_name, variables = c()) {
  
  crs <- get_crs(geom_data)
  
  geom_data <- st_coordinates(geom_data)
  
  hole_dif <- diff(geom_data[, 3])
  hole_dif[hole_dif < 0] <- 0
  part_dif <- diff(geom_data[, 4])
  part_dif[part_dif < 0] <- 0
  geom_dif <- diff(geom_data[, 5])
  
  geom_data[, 4] <- c(1, cumsum(hole_dif + part_dif + geom_dif) + 1)
  
  n_geoms <- length(unique(geom_data[, 5]))
  n_parts <- length(unique(geom_data[, 4]))
  
  nCoords <- nrow(geom_data)
  
  part_node_count <- rep(NA, n_parts)  
  part_type <- rep(pkg.env$multi_val, n_parts)
  node_count <- rep(NA, n_geoms)
  
  x_vals <- rep(NA, nCoords)
  y_vals <- x_vals
  nc_part <- 0
  coord <- 1
  
    for(geom in 1:n_geoms) {
      nd_count <- 0
      g_data <- geom_data[geom_data[, 5] == geom, ]
      
      nc_parts <- length(unique(g_data[, 4]))
      
      for(g_part in 1:nc_parts) {
        nc_part <- nc_part + 1
        
        if(g_part > 1 && g_data[g_data[, 4] == nc_part, 3][1] > 1) {
          part_type[nc_part] <- pkg.env$hole_val
        } else {
          part_type[nc_part] <-pkg.env$multi_val
        }
        
        coords <- g_data[g_data[, 4] == nc_part, c(1,2)]
        nrc <- nrow(coords)

        x_vals[coord:(coord + nrc - 1)] <- coords[1:nrc, 1]
        y_vals[coord:(coord + nrc - 1)] <- coords[1:nrc, 2]
        
        nd_count <- nd_count + nrc
        part_node_count[nc_part] <-  nrc
        coord <- coord + nrc
      }
      node_count[geom] <- nd_count
    }
  
  nc <- init_ncdf(nc_file, x_vals, y_vals, pkg.env$node_dim_name)
  
  write_geom_POLYGON(nc, node_count, instance_dim_name, part_node_count, part_type)
  
  put_meta_nc(nc, crs, variables)
  
  close.nc(nc)
  
  return(nc_file)
}

get_crs <- function(geom_data) {
  crs <- ncmeta::nc_prj_to_gridmapping(st_crs(geom_data)$proj4string)
  crs <- setNames(crs$value, crs$name)
  
  if(length(crs) == 0) {
    crs <- list(grid_mapping_name = "latitude_longitude",
                semi_major_axis = 6378137,
                inverse_flattening = 298.257223563,
                longitude_of_prime_meridian = 0)
    warning("No CRS was found. Assuming WGS84 Lat Lon.")
  }
  return(crs)
}

get_type <- function(geom_data) {
  type <- unique(st_geometry_type(geom_data))
  
  if(length(type) > 1) stop("Found multiple geometry types, only one is supported.")
  
  return(type)
}

init_ncdf <- function(nc_file, x_vals, y_vals, node_dim_name) {
  if(file.exists(nc_file)) {
    nc <- open.nc(nc_file, write = TRUE)
  } else {
    nc <- create.nc(nc_file, format = "classic4")
  }
  
  try(dim.def.nc(nc, node_dim_name, length(x_vals), unlim = FALSE), silent = TRUE)
  
  var.def.nc(nc, pkg.env$x_nodes, "NC_DOUBLE", node_dim_name)
  var.def.nc(nc, pkg.env$y_nodes, "NC_DOUBLE", node_dim_name)
  
  att.put.nc(nc, pkg.env$x_nodes, "units", "NC_CHAR", "degrees_east")
  att.put.nc(nc, pkg.env$y_nodes, "units", "NC_CHAR", "degrees_north")
  
  close.nc(nc)
  nc <- open.nc(nc_file, write = TRUE)
  
  var.put.nc(nc, pkg.env$x_nodes, x_vals)
  var.put.nc(nc, pkg.env$y_nodes, y_vals)
  # att.put.nc(nc, pkg.env$x_nodes, 'standard_name', "NC_CHAR", 'longitude')
  # att.put.nc(nc, pkg.env$y_nodes, 'standard_name', "NC_CHAR", 'latitude')
  att.put.nc(nc, pkg.env$x_nodes, "axis", "NC_CHAR", "X")
  att.put.nc(nc, pkg.env$y_nodes, "axis", "NC_CHAR", "Y")
  
  geom_container <- var.def.nc(nc, pkg.env$geom_container_var_name, "NC_INT", NA)
  att.put.nc(nc, pkg.env$geom_container_var_name, pkg.env$node_coordinates, "NC_CHAR", paste(pkg.env$x_nodes, pkg.env$y_nodes))
  
  return(nc)
}

put_meta_nc <- function(nc, crs, variables) {
  if(length(crs) > 0) {
    att.put.nc(nc, pkg.env$geom_container_var_name, pkg.env$crs, "NC_CHAR", pkg.env$crs_var_name)
    
    for(var in variables) {
      att.put.nc(nc, var, pkg.env$crs, "NC_CHAR", pkg.env$crs_var_name)
    }
    
    crs_var <- var.def.nc(nc, pkg.env$crs_var_name, vartype = "NC_INT", dimensions = NA)
    
    # close.nc(nc)
    # 
    # nc <- open.nc(nc_file, write = TRUE)
    
    for(crs_att in names(crs)) att.put.nc(nc, pkg.env$crs_var_name, crs_att, pkg.env$nc_types[class(crs[crs_att][[1]])][[1]], crs[crs_att][[1]])
  }
  
  for(var in variables) {
    att.put.nc(nc, var, pkg.env$geometry_container_att_name, "NC_CHAR", pkg.env$geom_container_var_name)
  }
  
  att.put.nc(nc, "NC_GLOBAL", "Conventions", "NC_CHAR", pkg.env$cf_version)
  
  return(invisible(nc))
}

write_geom_POINT <- function(nc) {
  att.put.nc(nc, pkg.env$geom_container_var_name, pkg.env$geom_type_attr_name, "NC_CHAR", 'point')
  return(invisible(nc))
}

write_geom_LINESTRING <- function(nc, node_count, instance_dim_name, part_node_count) {
  put_line(nc, node_count, instance_dim_name, part_node_count)
}

put_line <- function(nc, node_count, instance_dim_name, part_node_count) {
  att.put.nc(nc, pkg.env$geom_container_var_name, pkg.env$geom_type_attr_name, "NC_CHAR", 'line')
  
  put_node_count(nc, node_count, instance_dim_name)
  
  if(!is.null(part_node_count)) put_part_node_count(nc, part_node_count)
}

write_geom_POLYGON <- function(nc, node_count, instance_dim_name, part_node_count, part_type) {
  put_poly(nc, node_count, instance_dim_name, part_node_count, part_type)
}

put_poly <- function(nc, node_count, instance_dim_name, part_node_count, part_type) {
  att.put.nc(nc, pkg.env$geom_container_var_name, pkg.env$geom_type_attr_name, "NC_CHAR", 'polygon')
  
  put_node_count(nc, node_count, instance_dim_name)
  
  if(!is.null(part_node_count)) put_part_node_count(nc, part_node_count)
  
  if(!is.null(part_type)) put_hole_node_count(nc, part_type)
}

put_node_count <- function(nc, node_count, instance_dim_name) {
  att.put.nc(nc, pkg.env$geom_container_var_name, pkg.env$node_count_attr_name, "NC_CHAR", pkg.env$node_count_var_name)
  
  instance_dim <- tryCatch({
    dim.inq.nc(nc, instance_dim_name)},
    error = function(e) {
      dim.def.nc(nc, instance_dim_name, length(node_count), unlim = FALSE)
      dim.inq.nc(nc, instance_dim_name)
    })
  
  var.def.nc(nc, pkg.env$node_count_var_name, "NC_INT", instance_dim$id)
  att.put.nc(nc, pkg.env$node_count_var_name, "long_name", "NC_CHAR", "count of coordinates in each instance geometry")
  var.put.nc(nc, pkg.env$node_count_var_name, node_count)
  return(invisible(nc))
}

put_part_node_count <- function(nc, part_node_count) {
  dim.def.nc(nc, pkg.env$part_dim_name, length(part_node_count))
  var.def.nc(nc, pkg.env$part_node_count_var_name, "NC_INT", pkg.env$part_dim_name)
  var.put.nc(nc, pkg.env$part_node_count_var_name, part_node_count)
  att.put.nc(nc, pkg.env$part_node_count_var_name, "long_name", "NC_CHAR", "count of nodes in each geometry part")
  att.put.nc(nc, pkg.env$geom_container_var_name, pkg.env$part_node_count_attr_name, "NC_CHAR", pkg.env$part_node_count_var_name)
  return(invisible(nc))
}

put_hole_node_count <- function(nc, part_type) {
  var.def.nc(nc, pkg.env$part_type_var_name, "NC_INT", pkg.env$part_dim_name)
  var.put.nc(nc, pkg.env$part_type_var_name, part_type)
  att.put.nc(nc, pkg.env$part_type_var_name, "long_name", "NC_CHAR", "type of each geometry part")
  att.put.nc(nc, pkg.env$geom_container_var_name, pkg.env$part_type_attr_name, "NC_CHAR", pkg.env$part_type_var_name)
  return(invisible(nc))
}
USGS-R/netcdf.dsg documentation built on Feb. 9, 2024, 6:21 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close