Merge pull request #30 from dblodgett-usgs/master

NCDFSG Functionality

DESCRIPTION

@@ -1,17 +1,23 @@
 Package: netcdf.dsg
 Type: Package
 Title: Discrete Sampling Geometry NetCDF
-Version: 0.3.0
-Date: 2016-11-26
+Version: 0.4.0
+Date: 2017-07-17
 Author: Luke Winslow, Dave Blodgett
 Maintainer: David Blodgett <dblodgett@usgs.gov>
 Description: Tools to create discrete sampling geometry NetCDF files from R.
 Imports:
     ncdf4,
+    sp,
+  	rgdal,
     methods
+Depends:
+  R (>= 3.0)
 Suggests: 
 		testthat,
 		geoknife,
-		sp
+    maptools,
+    rgeos
 License: CC0
+LazyData: TRUE
 RoxygenNote: 6.0.1

NAMESPACE

@@ -1,5 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
+export(FromNCDFSG)
+export(ToNCDFSG)
+export(addGeomData)
+export(checkNCDF)
+export(getGmFromPrj)
+export(getPrjFromNCDF)
 export(read_instance_data)
 export(read_timeseries_dsg)
 export(write_instance_data)
@@ -18,3 +24,17 @@ importFrom(ncdf4,ncvar_add)
 importFrom(ncdf4,ncvar_def)
 importFrom(ncdf4,ncvar_get)
 importFrom(ncdf4,ncvar_put)
+importFrom(rgdal,CRSargs)
+importFrom(rgdal,checkCRSArgs)
+importFrom(sp,CRS)
+importFrom(sp,Line)
+importFrom(sp,Lines)
+importFrom(sp,Polygon)
+importFrom(sp,Polygons)
+importFrom(sp,SpatialLines)
+importFrom(sp,SpatialLinesDataFrame)
+importFrom(sp,SpatialPoints)
+importFrom(sp,SpatialPointsDataFrame)
+importFrom(sp,SpatialPolygons)
+importFrom(sp,SpatialPolygonsDataFrame)
+importFrom(sp,polygons)

R/AAA.R

@@ -1,4 +1,42 @@
 
 .onAttach <- function(libname, pkgname) {
-	packageStartupMessage("This information is preliminary or provisional and is subject to revision. It is being provided to meet the need for timely best science. The information has not received final approval by the U.S. Geological Survey (USGS) and is provided on the condition that neither the USGS nor the U.S. Government shall be held liable for any damages resulting from the authorized or unauthorized use of the information.")
-}
+	packageStartupMessage("This information is preliminary or provisional 
+and is subject to revision. It is being provided 
+to meet the need for timely best science. The 
+information has not received final approval by the 
+U.S. Geological Survey (USGS) and is provided on the 
+condition that neither the USGS nor the U.S. Government 
+shall be held liable for any damages resulting from the 
+authorized or unauthorized use of the information.
+
+****Support Package****
+This package is a USGS-R Support package. 
+see: https://owi.usgs.gov/R/packages.html#support")
+}
+
+pkg.env <- new.env()
+pkg.env$multi_val <- 0
+pkg.env$hole_val <- 1
+
+pkg.env$instance_dim_name <- "instance"
+
+# Arbitrary dim and variable names assumed in code.
+pkg.env$node_dim_name <- "node"
+pkg.env$part_dim_name <- "part"
+pkg.env$part_node_count_var_name <- "part_node_count"
+pkg.env$part_type_var_name <- "interior_ring"
+pkg.env$node_count_var_name <- "node_count"
+pkg.env$geom_container_var_name <- "geometry_container"
+pkg.env$crs_var_name <- "grid_mapping"
+
+# Variables prescribed in the specification.
+pkg.env$cf_version <- "CF-1.8"
+pkg.env$x_cf_role <- "geometry_x_node"
+pkg.env$y_cf_role <- "geometry_y_node"
+pkg.env$node_coordinates <- "node_coordinates"
+pkg.env$geom_type_attr_name <- "geometry_type"
+pkg.env$node_count_attr_name <- "node_count"
+pkg.env$part_node_count_attr_name <- "part_node_count"
+pkg.env$part_type_attr_name <- "interior_ring"
+pkg.env$geometry_container_att_name <- "geometry"
+pkg.env$crs <- "grid_mapping"

R/FromNCDFSG.R

@@ -0,0 +1,125 @@
+#'@title Convert NetCDF to sp objects
+#'
+#'
+#'@param nc_file A string file path to the nc file to be read.
+#'
+#'@description
+#'Attemps to convert a NetCDF-CF DSG Simple Geometry file into an sp object.
+#'
+#'@references
+#'https://github.com/twhiteaker/netCDF-CF-simple-geometry
+#'
+#'@importFrom ncdf4 nc_open nc_close ncvar_get ncatt_get
+#'@importFrom sp Polygon Polygons SpatialPolygons SpatialPolygonsDataFrame CRS Line Lines SpatialLines SpatialLinesDataFrame SpatialPointsDataFrame
+#'
+#'@return sp object containing spatial geometry of type found in the NetCDF-CF DSG file.
+#'
+#'@export
+FromNCDFSG = function(nc_file) {
+
+  nc <- nc_open(nc_file)
+
+  checkVals <- checkNCDF(nc)
+
+  instance_id<-checkVals$instance_id
+  instance_dim<-checkVals$instance_dim
+  geom_container <- checkVals$geom_container
+  variable_list <- checkVals$variable_list
+  crs <- checkVals$crs
+
+  line<-FALSE 
+  poly<-FALSE 
+  point<-FALSE
+  if(grepl("polygon", geom_container$geom_type)) { poly<-TRUE
+  } else if(grepl("line", geom_container$geom_type)) { line<-TRUE
+  } else point <- TRUE
+
+  xCoords <- c(ncvar_get(nc, geom_container$x))
+  yCoords <- c(ncvar_get(nc, geom_container$y))
+
+  if(length(crs) == 0) {
+    prj <- "+proj=longlat +datum=WGS84"
+  } else {
+    prj <- getPrjFromNCDF(crs)
+  }
+
+  if(point) {
+    point_data <- matrix(c(xCoords,
+                           yCoords), ncol=2)
+    dataFrame <- read_instance_data(nc, instance_dim)
+    if(geom_container$geom_type == "multipoint") {
+      stop("reading multipoint is not supported yet.")
+      # This is where handling for multipoint would go.
+    }
+    SPGeom <- SpatialPointsDataFrame(point_data, proj4string = CRS(prj),
+                                     data = dataFrame, match.ID = FALSE)
+  } else {
+    node_count <- c(ncvar_get(nc, geom_container$node_count))
+    if(!is.null(instance_id)) {
+      instance_names <- ncvar_get(nc, instance_id)
+    } else {
+      instance_names <- as.character(c(1:length(node_count)))
+    }
+    if(is.character(geom_container$part_node_count)) {
+      part_node_count <- ncvar_get(nc, geom_container$part_node_count)
+    } else {
+      part_node_count <- node_count
+    }
+    if(is.character(geom_container$part_type)) {
+      part_type <- ncvar_get(nc, geom_container$part_type)
+    } else {
+      part_type <- rep(pkg.env$multi_val, length(part_node_count))
+    }
+
+    node_start <- 1
+    geom_node_stop <- 0
+    pInd <- 1
+    Srl <- list()
+    for(geom in 1:length(node_count)) {
+
+      geom_node_stop <- geom_node_stop + node_count[geom]
+
+      srl <- list()
+
+      while(node_start < geom_node_stop) {
+        part_node_stop <- node_start + part_node_count[pInd] - 1
+
+        if(part_type[pInd] == pkg.env$hole_val) { hole <- TRUE
+        } else { hole <- FALSE }
+
+        coords <- matrix(c(xCoords[node_start:part_node_stop],yCoords[node_start:part_node_stop]),ncol=2)
+
+        if(poly) { tsrl<-Polygon(coords, hole=hole)
+        } else if(line) { tsrl<-Line(coords) }
+
+        dimnames(tsrl@coords) <- list(NULL, c("x", "y"))
+
+        srl <- append(srl, tsrl)
+
+        node_start <- node_start + part_node_count[pInd]; pInd <- pInd + 1
+      }
+      if(poly) {
+        Srl <- append(Srl, Polygons(srl, as.character(geom)))
+      }  else if(line) {
+        Srl <- append(Srl, Lines(srl, as.character(geom)))
+      }
+    }
+    dataFrame <- read_instance_data(nc, instance_dim)
+
+    for(varName in names(dataFrame)) {
+      if(!varName %in% variable_list) {
+        dataFrame[varName] <- NULL
+      }
+    }
+
+    if(poly) {
+      SPGeom <- SpatialPolygonsDataFrame(SpatialPolygons(Srl, proj4string = CRS(prj)),
+                                         dataFrame, match.ID = FALSE)
+    } else if(line) {
+      SPGeom <- SpatialLinesDataFrame(SpatialLines(Srl, proj4string = CRS(prj)),
+                                      dataFrame, match.ID = FALSE)
+    }
+  }
+  nc_close(nc)
+  return(SPGeom)
+}

R/ToNCDFSG.R

@@ -0,0 +1,84 @@
+#'@title Convert sp objects to NetCDF
+#'
+#'
+#'@param nc_file A string file path to the nc file to be created.
+#'@param instance_names A character vector of names for geometries.
+#'If NULL, integers are used. If the geomData has a data frame, this is not used.
+#'@param instance_dim_name If the file provided already has an instance dimension,
+#'it needs to be provided as a character string otherwise a new instance dim may be created.
+#'@param geomData An object of class \code{SpatialPoints}, \code{SpatialLines} or
+#'\code{SpatialPolygons} with WGS84 lon in the x coordinate and lat in the y coordinate.
+#'Note that three dimensional geometries is not supported.
+#'@param lats Vector of WGS84 latitudes
+#'@param lons Vector of WGS84 longitudes
+#'@param variables If a an existing netcdf files is provided, this list of strings is used
+#'to add the geometry container attribute to the named existing variables.
+#'
+#'@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 ncdf4 nc_open ncvar_add nc_close ncvar_def ncvar_put ncatt_put ncdim_def
+#'@importFrom sp SpatialLinesDataFrame polygons SpatialPoints
+#'
+#'@export
+ToNCDFSG = function(nc_file, geomData = NULL, instance_names = NULL, instance_dim_name = NULL, lats = NULL, lons = NULL, variables = list()){
+
+  pointsMode <- FALSE
+
+  if(is.null(instance_names) && !is.null(geomData)) {
+    if(class(geomData)=="SpatialPoints" || class(geomData)=="SpatialPointsDataFrame") {
+      instance_names <- as.character(unique(attributes(geomData@coords)$dimnames[[1]]))
+    } else {
+      instance_names <- as.character(c(1:length(geomData)))
+    }
+  }
+
+  if(class(geomData) == "SpatialPolygonsDataFrame") {
+    attData<-geomData@data
+    geomData<-polygons(geomData)
+  } else if(class(geomData) == "SpatialLinesDataFrame") {
+    attData<-geomData@data
+  } else if(class(geomData) == "SpatialPolygons") {
+    geomData<-polygons(geomData)
+  } else if(class(geomData) == "SpatialLines") {
+    geomData<-SpatialLinesDataFrame(geomData,data=as.data.frame(instance_names,stringsAsFactors = FALSE))
+  } else if(class(geomData) == "SpatialPoints") {
+    pointsMode<-TRUE
+  } else if(class(geomData) == "SpatialPointsDataFrame") {
+    pointsMode<-TRUE
+    attData<-geomData@data
+  } else if(!is.null(lats)) {
+    pointsMode<-TRUE
+    geomData <- SpatialPoints(as.data.frame(list(x=lons, y=lats)),proj4string = CRS("+proj=longlat +datum=WGS84"))
+    if(is.null(instance_names)) {
+      instance_names<-as.character(c(1:length(lats)))
+    }
+  } else {
+    stop("Did not find supported spatial data.")
+  }
+
+  if(!pointsMode && !is.null(geomData)) {
+    if(length(instance_names)!=length(geomData)) stop('instance_names must be same length as data')
+  }
+
+  if(is.null(instance_dim_name)) {
+    instance_dim_name <- pkg.env$instance_dim_name
+  }
+
+  if(exists("attData")) {
+    itemp <- sapply(attData, is.factor)
+    attData[itemp] <- lapply(attData[itemp], as.character)
+    nc_file <- write_instance_data(nc_file, attData, instance_dim_name)
+    variables <- c(variables, names(attData))
+  }
+
+  nc_file <- addGeomData(nc_file, geomData, instance_dim_name, variables = variables)
+
+  return(nc_file)
+}
+
+