Reading and Writing NetCDF-CF Geometry

Source: vignettes/geometry.Rmd

Introduction

This example shows how to

  1. add polygon data to a NetCDF-CF Discrete Sampling Geometry timeSeries featureType file,
  2. read the polygon data from the NetCDF file,
  3. and write it to a standard GIS format.

Below the example, all the geometry types that ncdfgeom handles are shown in detail.

Load Spatial Data

Tables of point, line, or polygon features with associated timeseries are the target for this functionality. Here, we use some sample data from the ncdfgeom package.

example_file <- tempfile()

file.copy(from = system.file('extdata/example_huc_eta.nc', package = 'ncdfgeom'), 
                    to = example_file, 
                    overwrite = TRUE) -> quiet

polygons <- sf::read_sf(system.file('extdata/example_huc_eta.json', package = 'ncdfgeom'))

polygons <- dplyr::select(polygons, LOADDATE, AREASQKM, HUC12, NAME)

plot(sf::st_geometry(polygons))

Now we have the polygons as shown above and a NetCDF file with a header that looks like:

Now we can use the write_geometry function to add the polygon data to the NetCDF file.

(vars <- ncmeta::nc_vars(example_file))
#> # A tibble: 5 × 5
#>      id name         type      ndims natts
#>   <int> <chr>        <chr>     <int> <int>
#> 1     0 lat          NC_DOUBLE     1     4
#> 2     1 lon          NC_DOUBLE     1     4
#> 3     2 time         NC_DOUBLE     1     4
#> 4     3 station_name NC_CHAR       2     5
#> 5     4 et           NC_INT        2     4

ncdfgeom::write_geometry(nc_file=example_file,
                         geom_data = polygons, 
                         instance_dim_name = "station", 
                         variables = vars$name) -> example_file

Now the NetCDF file looks like:

#> netcdf file63806c813359 {
#> dimensions:
#>  maxStrlen64 = 64 ;
#>  station = 2 ;
#>  time = 25 ;
#>  char = 12 ;
#>  node = 2628 ;
#> variables:
#>  double lat(station) ;
#>      lat:units = "degrees_north" ;
#>      lat:missing_value = -999. ;
#>      lat:long_name = "latitude of the observation" ;
#>      lat:standard_name = "latitude" ;
#>      lat:grid_mapping = "grid_mapping" ;
#>      lat:geometry = "geometry_container" ;
#>  double lon(station) ;
#>      lon:units = "degrees_east" ;
#>      lon:missing_value = -999. ;
#>      lon:long_name = "longitude of the observation" ;
#>      lon:standard_name = "longitude" ;
#>      lon:grid_mapping = "grid_mapping" ;
#>      lon:geometry = "geometry_container" ;
#>  double time(time) ;
#>      time:units = "days since 1970-01-01 00:00:00" ;
#>      time:missing_value = -999. ;
#>      time:long_name = "time of measurement" ;
#>      time:standard_name = "time" ;
#>      time:grid_mapping = "grid_mapping" ;
#>      time:geometry = "geometry_container" ;
#>  char station_name(station, maxStrlen64) ;
#>      station_name:units = "" ;
#>      station_name:missing_value = "" ;
#>      station_name:long_name = "Station Names" ;
#>      station_name:cf_role = "timeseries_id" ;
#>      station_name:standard_name = "station_id" ;
#>      station_name:grid_mapping = "grid_mapping" ;
#>      station_name:geometry = "geometry_container" ;
#>  int et(station, time) ;
#>      et:units = "mm" ;
#>      et:missing_value = -999 ;
#>      et:long_name = "Area Weighted Mean Actual Evapotranspiration" ;
#>      et:coordinates = "time lat lon" ;
#>      et:grid_mapping = "grid_mapping" ;
#>      et:geometry = "geometry_container" ;
#>  char LOADDATE(station, char) ;
#>      LOADDATE:units = "unknown" ;
#>      LOADDATE:grid_mapping = "grid_mapping" ;
#>      LOADDATE:geometry = "geometry_container" ;
#>  double AREASQKM(station) ;
#>      AREASQKM:units = "unknown" ;
#>      AREASQKM:missing_value = -9999.999 ;
#>      AREASQKM:grid_mapping = "grid_mapping" ;
#>      AREASQKM:geometry = "geometry_container" ;
#>  char HUC12(station, char) ;
#>      HUC12:units = "unknown" ;
#>      HUC12:grid_mapping = "grid_mapping" ;
#>      HUC12:geometry = "geometry_container" ;
#>  char NAME(station, char) ;
#>      NAME:units = "unknown" ;
#>      NAME:grid_mapping = "grid_mapping" ;
#>      NAME:geometry = "geometry_container" ;
#>  double x_nodes(node) ;
#>      x_nodes:units = "degrees_east" ;
#>      x_nodes:axis = "X" ;
#>  double y_nodes(node) ;
#>      y_nodes:units = "degrees_north" ;
#>      y_nodes:axis = "Y" ;
#>  int geometry_container ;
#>      geometry_container:node_coordinates = "x_nodes y_nodes" ;
#>      geometry_container:geometry_type = "polygon" ;
#>      geometry_container:node_count = "node_count" ;
#>      geometry_container:grid_mapping = "grid_mapping" ;
#>  int node_count(station) ;
#>      node_count:long_name = "count of coordinates in each instance geometry" ;
#>  int grid_mapping ;
#>      grid_mapping:grid_mapping_name = "latitude_longitude" ;
#>      grid_mapping:semi_major_axis = 6378137. ;
#>      grid_mapping:inverse_flattening = 298.257223563 ;
#>      grid_mapping:longitude_of_prime_meridian = 0. ;
#> 
#> // global attributes:
#>      :Conventions = "CF-1.8" ;
#>      :featureType = "timeSeries" ;
#>      :cdm_data_type = "Station" ;
#> }

Read the polygon data from the file and write it out to a geopackage.

polygons_sf <- ncdfgeom::read_geometry(example_file)

plot(sf::st_geometry(polygons_sf))

sf::write_sf(polygons_sf, "polygons.gpkg")

Geometry Types

Examples - Contiguous Ragged Arrays

Point (2D)

Well Known Text (WKT): POINT (30 10)
Common Data Language (CDL):

netcdf sample_point {
dimensions:
    instance = 1 ;
variables:
    double x_nodes(instance) ;
        x_nodes:units = "degrees_east" ;
        x_nodes:axis = "X" ;
    double y_nodes(instance) ;
        y_nodes:units = "degrees_north" ;
        y_nodes:axis = "Y" ;
    int geometry_container ;
        geometry_container:node_coordinates = "x_nodes y_nodes" ;
        geometry_container:geometry_type = "point" ;
        geometry_container:grid_mapping = "grid_mapping" ;
    int grid_mapping ;
        grid_mapping:grid_mapping_name = "latitude_longitude" ;
        grid_mapping:semi_major_axis = 6378137. ;
        grid_mapping:inverse_flattening = 298.257223563 ;
        grid_mapping:longitude_of_prime_meridian = 0. ;

// global attributes:
        :Conventions = "CF-1.8" ;
data:

 x_nodes = 30 ;

 y_nodes = 10 ;

 geometry_container = _ ;

 grid_mapping = _ ;
}

LineString (2D)

Well Known Text (WKT): LINESTRING (30 10, 10 30, 40 40)
Common Data Language (CDL):

netcdf sample_linestring {
dimensions:
    node = 3 ;
    instance = 1 ;
variables:
    double x_nodes(node) ;
        x_nodes:units = "degrees_east" ;
        x_nodes:axis = "X" ;
    double y_nodes(node) ;
        y_nodes:units = "degrees_north" ;
        y_nodes:axis = "Y" ;
    int geometry_container ;
        geometry_container:node_coordinates = "x_nodes y_nodes" ;
        geometry_container:geometry_type = "line" ;
        geometry_container:node_count = "node_count" ;
        geometry_container:grid_mapping = "grid_mapping" ;
    int node_count(instance) ;
        node_count:long_name = "count of coordinates in each instance geometry" ;
    int grid_mapping ;
        grid_mapping:grid_mapping_name = "latitude_longitude" ;
        grid_mapping:semi_major_axis = 6378137. ;
        grid_mapping:inverse_flattening = 298.257223563 ;
        grid_mapping:longitude_of_prime_meridian = 0. ;

// global attributes:
        :Conventions = "CF-1.8" ;
data:

 x_nodes = 30, 10, 40 ;

 y_nodes = 10, 30, 40 ;

 geometry_container = _ ;

 node_count = 3 ;

 grid_mapping = _ ;
}

Polygon (2D)

Well Known Text (WKT): POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))
Common Data Language (CDL):

netcdf sample_polygon {
dimensions:
    node = 5 ;
    instance = 1 ;
variables:
    double x_nodes(node) ;
        x_nodes:units = "degrees_east" ;
        x_nodes:axis = "X" ;
    double y_nodes(node) ;
        y_nodes:units = "degrees_north" ;
        y_nodes:axis = "Y" ;
    int geometry_container ;
        geometry_container:node_coordinates = "x_nodes y_nodes" ;
        geometry_container:geometry_type = "polygon" ;
        geometry_container:node_count = "node_count" ;
        geometry_container:grid_mapping = "grid_mapping" ;
    int node_count(instance) ;
        node_count:long_name = "count of coordinates in each instance geometry" ;
    int grid_mapping ;
        grid_mapping:grid_mapping_name = "latitude_longitude" ;
        grid_mapping:semi_major_axis = 6378137. ;
        grid_mapping:inverse_flattening = 298.257223563 ;
        grid_mapping:longitude_of_prime_meridian = 0. ;

// global attributes:
        :Conventions = "CF-1.8" ;
data:

 x_nodes = 30, 40, 20, 10, 30 ;

 y_nodes = 10, 40, 40, 20, 10 ;

 geometry_container = _ ;

 node_count = 5 ;

 grid_mapping = _ ;
}

MultiLineString (2D)

Well Known Text (WKT): MULTILINESTRING ((10 10, 20 20, 10 40), (40 40, 30 30, 40 20, 30 10))
Common Data Language (CDL):

netcdf sample_multilinestring {
dimensions:
    node = 7 ;
    instance = 1 ;
    part = 2 ;
variables:
    double x_nodes(node) ;
        x_nodes:units = "degrees_east" ;
        x_nodes:axis = "X" ;
    double y_nodes(node) ;
        y_nodes:units = "degrees_north" ;
        y_nodes:axis = "Y" ;
    int geometry_container ;
        geometry_container:node_coordinates = "x_nodes y_nodes" ;
        geometry_container:geometry_type = "line" ;
        geometry_container:node_count = "node_count" ;
        geometry_container:part_node_count = "part_node_count" ;
        geometry_container:grid_mapping = "grid_mapping" ;
    int node_count(instance) ;
        node_count:long_name = "count of coordinates in each instance geometry" ;
    int part_node_count(part) ;
        part_node_count:long_name = "count of nodes in each geometry part" ;
    int grid_mapping ;
        grid_mapping:grid_mapping_name = "latitude_longitude" ;
        grid_mapping:semi_major_axis = 6378137. ;
        grid_mapping:inverse_flattening = 298.257223563 ;
        grid_mapping:longitude_of_prime_meridian = 0. ;

// global attributes:
        :Conventions = "CF-1.8" ;
data:

 x_nodes = 10, 20, 10, 40, 30, 40, 30 ;

 y_nodes = 10, 20, 40, 40, 30, 20, 10 ;

 geometry_container = _ ;

 node_count = 7 ;

 part_node_count = 3, 4 ;

 grid_mapping = _ ;
}

MultiPolygon (2D)

Well Known Text (WKT): MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)), ((15 5, 40 10, 10 20, 5 10, 15 5)))
Common Data Language (CDL):

netcdf sample_multipolygon {
dimensions:
    node = 9 ;
    instance = 1 ;
    part = 2 ;
variables:
    double x_nodes(node) ;
        x_nodes:units = "degrees_east" ;
        x_nodes:axis = "X" ;
    double y_nodes(node) ;
        y_nodes:units = "degrees_north" ;
        y_nodes:axis = "Y" ;
    int geometry_container ;
        geometry_container:node_coordinates = "x_nodes y_nodes" ;
        geometry_container:geometry_type = "polygon" ;
        geometry_container:node_count = "node_count" ;
        geometry_container:interior_ring = "interior_ring" ;
        geometry_container:grid_mapping = "grid_mapping" ;
        geometry_container:part_node_count = "part_node_count" ;
    int node_count(instance) ;
        node_count:long_name = "count of coordinates in each instance geometry" ;
    int part_node_count(part) ;
        part_node_count:long_name = "count of nodes in each geometry part" ;
    int interior_ring(part) ;
        interior_ring:long_name = "type of each geometry part" ;
    int grid_mapping ;
        grid_mapping:grid_mapping_name = "latitude_longitude" ;
        grid_mapping:semi_major_axis = 6378137. ;
        grid_mapping:inverse_flattening = 298.257223563 ;
        grid_mapping:longitude_of_prime_meridian = 0. ;

// global attributes:
        :Conventions = "CF-1.8" ;
data:

 x_nodes = 30, 45, 10, 30, 15, 40, 10, 5, 15 ;

 y_nodes = 20, 40, 40, 20, 5, 10, 20, 10, 5 ;

 geometry_container = _ ;

 node_count = 9 ;

 part_node_count = 4, 5 ;

 interior_ring = 0, 0 ;

 grid_mapping = _ ;
}

Polygon with One Interior Ring (2D)

Well Known Text (WKT): POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10), (20 30, 35 35, 30 20, 20 30))
Common Data Language (CDL):

netcdf sample_polygon_hole {
dimensions:
    node = 9 ;
    instance = 1 ;
    part = 2 ;
variables:
    double x_nodes(node) ;
        x_nodes:units = "degrees_east" ;
        x_nodes:axis = "X" ;
    double y_nodes(node) ;
        y_nodes:units = "degrees_north" ;
        y_nodes:axis = "Y" ;
    int geometry_container ;
        geometry_container:node_coordinates = "x_nodes y_nodes" ;
        geometry_container:geometry_type = "polygon" ;
        geometry_container:node_count = "node_count" ;
        geometry_container:interior_ring = "interior_ring" ;
        geometry_container:grid_mapping = "grid_mapping" ;
        geometry_container:part_node_count = "part_node_count" ;
    int node_count(instance) ;
        node_count:long_name = "count of coordinates in each instance geometry" ;
    int part_node_count(part) ;
        part_node_count:long_name = "count of nodes in each geometry part" ;
    int interior_ring(part) ;
        interior_ring:long_name = "type of each geometry part" ;
    int grid_mapping ;
        grid_mapping:grid_mapping_name = "latitude_longitude" ;
        grid_mapping:semi_major_axis = 6378137. ;
        grid_mapping:inverse_flattening = 298.257223563 ;
        grid_mapping:longitude_of_prime_meridian = 0. ;

// global attributes:
        :Conventions = "CF-1.8" ;
data:

 x_nodes = 35, 45, 15, 10, 35, 20, 35, 30, 20 ;

 y_nodes = 10, 45, 40, 20, 10, 30, 35, 20, 30 ;

 geometry_container = _ ;

 node_count = 9 ;

 part_node_count = 5, 4 ;

 interior_ring = 0, 1 ;

 grid_mapping = _ ;
}

MultiPolygon with One Interior Ring (2D)

Well Known Text (WKT): MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)), ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20)))
Common Data Language (CDL):

netcdf sample_multipolygon_hole {
dimensions:
    node = 14 ;
    instance = 1 ;
    part = 3 ;
variables:
    double x_nodes(node) ;
        x_nodes:units = "degrees_east" ;
        x_nodes:axis = "X" ;
    double y_nodes(node) ;
        y_nodes:units = "degrees_north" ;
        y_nodes:axis = "Y" ;
    int geometry_container ;
        geometry_container:node_coordinates = "x_nodes y_nodes" ;
        geometry_container:geometry_type = "polygon" ;
        geometry_container:node_count = "node_count" ;
        geometry_container:interior_ring = "interior_ring" ;
        geometry_container:grid_mapping = "grid_mapping" ;
        geometry_container:part_node_count = "part_node_count" ;
    int node_count(instance) ;
        node_count:long_name = "count of coordinates in each instance geometry" ;
    int part_node_count(part) ;
        part_node_count:long_name = "count of nodes in each geometry part" ;
    int interior_ring(part) ;
        interior_ring:long_name = "type of each geometry part" ;
    int grid_mapping ;
        grid_mapping:grid_mapping_name = "latitude_longitude" ;
        grid_mapping:semi_major_axis = 6378137. ;
        grid_mapping:inverse_flattening = 298.257223563 ;
        grid_mapping:longitude_of_prime_meridian = 0. ;

// global attributes:
        :Conventions = "CF-1.8" ;
data:

 x_nodes = 40, 20, 45, 40, 20, 10, 10, 30, 45, 20, 30, 20, 20, 30 ;

 y_nodes = 40, 45, 30, 40, 35, 30, 10, 5, 20, 35, 20, 15, 25, 20 ;

 geometry_container = _ ;

 node_count = 14 ;

 part_node_count = 4, 6, 4 ;

 interior_ring = 0, 0, 1 ;

 grid_mapping = _ ;
}

Multiple MultiPolygons with Interior Rings (2D)

Well Known Text (WKT): MULTIPOLYGON(((0 0, 20 0, 20 20, 0 20, 0 0), (1 1, 10 5, 19 1, 1 1), (5 15, 7 19, 9 15, 5 15), (11 15, 13 19, 15 15, 11 15)), ((5 25, 9 25, 7 29, 5 25)), ((11 25, 15 25, 13 29, 11 25)))
Common Data Language (CDL):

netcdf sample_multipolygons_holes {
dimensions:
    node = 25 ;
    instance = 1 ;
    part = 6 ;
variables:
    double x_nodes(node) ;
        x_nodes:units = "degrees_east" ;
        x_nodes:axis = "X" ;
    double y_nodes(node) ;
        y_nodes:units = "degrees_north" ;
        y_nodes:axis = "Y" ;
    int geometry_container ;
        geometry_container:node_coordinates = "x_nodes y_nodes" ;
        geometry_container:geometry_type = "polygon" ;
        geometry_container:node_count = "node_count" ;
        geometry_container:interior_ring = "interior_ring" ;
        geometry_container:grid_mapping = "grid_mapping" ;
        geometry_container:part_node_count = "part_node_count" ;
    int node_count(instance) ;
        node_count:long_name = "count of coordinates in each instance geometry" ;
    int part_node_count(part) ;
        part_node_count:long_name = "count of nodes in each geometry part" ;
    int interior_ring(part) ;
        interior_ring:long_name = "type of each geometry part" ;
    int grid_mapping ;
        grid_mapping:grid_mapping_name = "latitude_longitude" ;
        grid_mapping:semi_major_axis = 6378137. ;
        grid_mapping:inverse_flattening = 298.257223563 ;
        grid_mapping:longitude_of_prime_meridian = 0. ;

// global attributes:
        :Conventions = "CF-1.8" ;
data:

 x_nodes = 0, 20, 20, 0, 0, 1, 10, 19, 1, 5, 7, 9, 5, 11, 13, 15, 11, 5, 9, 
    7, 5, 11, 15, 13, 11 ;

 y_nodes = 0, 0, 20, 20, 0, 1, 5, 1, 1, 15, 19, 15, 15, 15, 19, 15, 15, 25, 
    25, 29, 25, 25, 25, 29, 25 ;

 geometry_container = _ ;

 node_count = 25 ;

 part_node_count = 5, 4, 4, 4, 4, 4 ;

 interior_ring = 0, 1, 1, 1, 0, 0 ;

 grid_mapping = _ ;
}