Merge c1ede8a into 33b464b

DESCRIPTION

@@ -1,15 +1,15 @@
 Package: ncdfgeom
 Type: Package
 Title: 'NetCDF' Geometry and Time Series
-Version: 1.0.0
+Version: 1.1.0
 Date: 2019-06-05
 Authors@R: c(person("David", "Blodgett", role = c("aut", "cre"), 
                     email = "dblodgett@usgs.gov"),
              person("Luke", "Winslow", role = "ctb"))
 Description: Tools to create time series and geometry 'NetCDF' files.
 URL: https://code.usgs.gov/water/ncdfgeom
 BugReports: https://github.com/USGS-R/ncdfgeom/issues
-Imports: RNetCDF, ncmeta, sf, dplyr, methods
+Imports: RNetCDF, ncmeta, sf, dplyr, methods, stars
 Depends:
   R (>= 3.0)
 Suggests: testthat, knitr, rmarkdown, pkgdown, tidyverse, sp, geoknife, ncdf4, jsonlite

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(st_as_stars,ncdfgeom)
 export(read_attribute_data)
 export(read_geometry)
 export(read_timeseries_dsg)
@@ -37,4 +38,6 @@ importFrom(sf,st_polygon)
 importFrom(sf,st_set_geometry)
 importFrom(sf,st_sf)
 importFrom(sf,st_sfc)
+importFrom(sf,st_zm)
+importFrom(stars,st_as_stars)
 importFrom(stats,setNames)

R/read_timeseries_dsg.R

@@ -197,5 +197,8 @@ read_timeseries_dsg = function(nc_file){
 	                                                    name == "processing_level")$value
 	nc_list$global_attributes$nc_title <- filter(nc_atts, variable == "NC_GLOBAL" &
 	                                               name == "title")$value
+
+	attr(nc_list, "class") <- "ncdfgeom"
+
 	return(nc_list)
 }

R/st_as_stars.R

@@ -0,0 +1,43 @@
+#' Convert ncdfgeom object into stars object.
+#' @importFrom stars st_as_stars
+#' @param .x Object of class ncdfgeom as returned by read_timeseries_dsg.
+#' @param ... not used.
+#' @param sf_geometry sf data.frame with geometry and attributes to be added to stars object. 
+#' Must have same number of rows as timeseries instances.
+#' @name st_as_stars
+#' @export
+#' 
+st_as_stars.ncdfgeom <- function(.x, ..., sf_geometry = NA) {
+  crs <- st_crs(4326)$proj4string
+  ts_points <- data.frame(X = .x$lons, Y = .x$lats, Z = .x$alts)
+  ts_points <- sf::st_as_sf(ts_points, coords = c("X", "Y", "Z"), crs = crs)
+
+  data <- .x$data_frames[[1]]
+  # data[["T"]] <- .x$time
+
+  gdim <- stars:::create_dimension(from = 1, to = length(.x$lats), 
+                                   refsys = crs, point = TRUE, 
+                                   values = ts_points$geometry)
+  tdim <- stars:::create_dimension(from = 1, to = length(.x$time), 
+                                   refsys = "POSIXct", point = FALSE, 
+                                   values = as.POSIXct(.x$time))
+  dim <- list(time = tdim, points = gdim)
+
+  if("sf" %in% class(sf_geometry)) {
+    if(length(gdim$values) != length(st_geometry(sf_geometry))) 
+      stop("geometry must be same length as instance dimension of timeseries")
+
+    is_point <- any(grepl("point", class(st_geometry(sf_geometry)), ignore.case = TRUE))
+
+    sf_dim <- stars:::create_dimension(from = 1, to = length(gdim$values),
+                                       refsys = st_crs(sf_geometry)$proj4string, 
+                                       point = is_point, is_raster = FALSE,
+                                       values = st_geometry(sf_geometry))
+
+    dim <- c(dim, list(geometry = sf_dim))
+  }
+
+  stars:::st_stars(x = setNames(list(as.matrix(.x$data_frames[[1]])), 
+                                .x$varmeta[[1]]$name), 
+                   dimensions =  stars:::create_dimensions(dim))
+}

R/write_attribute_data.R

@@ -57,7 +57,7 @@ write_attribute_data <- function(nc_file, att_data, instance_dim_name = "instanc
 	types <- list(numeric="NC_DOUBLE", integer = "NC_INT", character="NC_CHAR")
 
 	# Convert any dates to character. This could be improved later.
-	i <- sapply(att_data, is, class2 = "Date")
+	i <- sapply(att_data, is, class2 = "Date") | sapply(att_data, is, class2 = "POSIXt")
 	att_data[i] <- lapply(att_data[i], as.character)
 
 	charDimLen<-0

R/write_geometry.R

@@ -89,9 +89,9 @@ write_geometry = function(nc_file, geom_data, instance_dim_name = NULL, variable
 #'
 #' @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
+#' @importFrom sf st_geometry_type st_crs st_coordinates st_zm
 #' @noRd
-write_geom_data <- function(geom_data, ...) 
+write_geom_data <- function(geom_data, ...)
   UseMethod("write_geom_data")
 
 #' @noRd
@@ -159,6 +159,11 @@ write_geom_data.sfc_LINESTRING <- function(geom_data, nc_file, instance_dim_name
 #' @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)
@@ -186,7 +191,7 @@ write_geom_data.sfc_MULTILINESTRING <- function(geom_data, nc_file,
 
     g_data <- geom_data[geom_data[, 4] == geom, ]
 
-    for(g_part in 1:length(unique(g_data[, 3]))) {
+    for(g_part in unique(g_data[, 3])) {
       nc_part <- nc_part + 1
 
       coords <- g_data[g_data[, 3] == g_part, c(1,2)]

tests/testthat/helper-functions.R

@@ -97,5 +97,36 @@ get_sample_timeseries_data <- function() {
               lons = lons,
               lats = lats,
               alts = alts,
-              units = units))
+              units = units,
+              geom = yahara))
+}
+
+get_test_ncdf_object <- function(nc_file = tempfile()) {
+  nc_summary<-'test summary'
+  nc_date_create<-'2099-01-01'
+  nc_creator_name='test creator'
+  nc_creator_email='test@test.com'
+  nc_project='testthat ncdfgeom'
+  nc_proc_level='just a test no processing'
+  nc_title<-'test title'
+  global_attributes<-list(title = nc_title, summary = nc_summary, date_created=nc_date_create, 
+                          creator_name=nc_creator_name,creator_email=nc_creator_email,
+                          project=nc_project, processing_level=nc_proc_level)
+
+  test_data <- get_sample_timeseries_data()
+
+  testnc<-write_timeseries_dsg(nc_file, 
+                               names(test_data$var_data), 
+                               test_data$lats, test_data$lons, 
+                               as.character(test_data$time), 
+                               test_data$var, 
+                               test_data$alts, 
+                               data_unit=test_data$units,	
+                               data_prec='double',
+                               data_metadata=test_data$meta,
+                               attributes=global_attributes)
+
+  test_nc <- write_geometry(nc_file, test_data$geom, variables = test_data$meta$name)
+
+  list(ncdfgeom = read_timeseries_dsg(nc_file), sf = read_geometry(nc_file))
 }

tests/testthat/test_line.R

@@ -64,3 +64,17 @@ test_that("shapefile line data works", {
   }
 })
 
+test_that("NHDPlus Multilinestring", {
+  f <- system.file("extdata/nhdp_flowline_sample.gpkg", package = "ncdfgeom")
+
+  test_dat <- sf::read_sf(f)
+
+  test_nc <- expect_warning(write_geometry(tempfile(), test_dat), "Found more than two dimensions in geometry. Removing Z and M content.")
+
+  test_dat_2 <- read_geometry(test_nc)
+
+  expect_equal(class(test_dat_2$FDATE), "character") # coerced to character
+  expect_equal(class(sf::st_geometry(test_dat_2)[[1]])[1], "XY") 
+
+})
+

tests/testthat/test_read-write_timeseries_dsg.R

@@ -2,6 +2,7 @@ context("orthogonal netcdf timeseries")
 
 test_that("Create basic DSG file", {
 
+  # NOTE: this code has been moved to helper files but was left here to not mess with it.
   nc_file<-tempfile()
   nc_summary<-'test summary'
   nc_date_create<-'2099-01-01'
@@ -222,6 +223,7 @@ test_that("Create basic DSG file", {
 	expect_equivalent(testlist$global_attributes$nc_proc_level,'just a test no processing')
 	expect_equivalent(testlist$global_attributes$nc_title,'test title')
 	expect_equivalent(testlist$data_frames[1][[1]],test_data$var_data) # Plan to have the dataframes work for 1 to many variables.
+	expect_s3_class(testlist, "ncdfgeom")
 
   nc_file_borked <- tempfile()
   file.copy(nc_file, nc_file_borked, overwrite = TRUE)

tests/testthat/test_st_as_stars.R

@@ -0,0 +1,27 @@
+context("st_as_stars tests")
+
+test_that("basic st_as_stars", {
+  test_list <- get_test_ncdf_object()
+
+  stars_obj <- st_as_stars(test_list$ncdfgeom)
+
+  expect_s3_class(stars_obj, "stars")
+
+  dim <- stars::st_dimensions(stars_obj)
+  expect_equal(sf::st_crs(dim$points$refsys), sf::st_crs(4326))
+  expect_equal(dim$time$refsys, "POSIXct")
+
+  expect_s3_class(dim$points$values, "sfc_POINT")
+
+  expect_true(dim$points$point)  
+
+  stars_obj <- st_as_stars(test_list$ncdfgeom, sf_geometry = test_list$sf)
+
+  dim <- stars::st_dimensions(stars_obj)
+  expect_equal(sf::st_crs(dim$geometry$refsys), sf::st_crs(test_list$sf))
+
+  expect_s3_class(dim$geometry$values, "sfc_POLYGON")
+
+  expect_false(dim$geometry$point)  
+
+})