1

I have NetCDF files (e.g https://data.ceda.ac.uk/neodc/esacci/lakes/data/lake_products/L3S/v1.0/2019 global domain), and I want to extract the data based on a shapefile boundary ( in this case a Lake here - https://www.sciencebase.gov/catalog/item/530f8a0ee4b0e7e46bd300dd) and then save clipped data as a NetCDF file but retain all the original metadata and variables names within the clipped file. This is what I have done far

library(rgdal)
library(sf)
library(ncdf4)
library(terra)

#Read in the shapefile of Lake 
Lake_shape <- readOGR("C:/Users/CEDA/hydro_p_LakeA/hydro_p_A.shp")
# Reading the netcdf file using Terra Package function rast
test <- rast("ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190705-fv1.0.nc")
# List of some of variables names for orginal dataset 
      head(names(test))
[1] "water_surface_height_above_reference_datum" "water_surface_height_uncertainty"           "lake_surface_water_extent"                 
[4] "lake_surface_water_extent_uncertainty"      "lake_surface_water_temperature"             "lswt_uncertainty"   
                                 
#Clipping data to smaller Lake domain using the crop function in Terra Package
test3 <- crop(test, Lake_shape)
#Listing the some variables names for clipped data
head(names(test3))
[1] "water_surface_height_above_reference_datum" "water_surface_height_uncertainty"           "lake_surface_water_extent"                 
[4] "lake_surface_water_extent_uncertainty"      "lake_surface_water_temperature"             "lswt_uncertainty" 


# Writing the crop dataset as netcdf or Raster Layer using the WriteCDF function 

filepath<-"Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0"
fname <- paste0( "C:/Users/CEDA/",filepath,".nc")
rnc <- writeCDF(test3, filename =fname, overwrite=T)”

My main issue here when I read in clipped the netCDF file I don’t seem to be able to keep the names of the data variables of the original NetCDF. They are all being renamed automatically when I am saving the clipped dataset as a new netCDF using the writeCDF function.

#Reading in the new clipped file
 LakeA<-rast("Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0.nc")
> head(names(LakeA))
[1] "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_1" "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_2"
[3] "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_3" "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_4"
[5] "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_5" "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_6"

So is it possible to clone/copy all the metadata variables from the original NetCDF dataset when clipping to the smaller domain/shapefile in R, then saving as NetCDF? Any guidance on how to do this in R would be really appreciated. (NetCDF and R are all new to me so I am not sure what I am missing or have the in-depth knowledge to sort this).

3
  • 1
    try: rnc <- writeCDF(test3, filename =fname,varname=names(test3), overwrite=T) So that the layer names are assigned by the names in the SpatRaster. It should output the names you are after when reading It back in
    – Stackbeans
    Commented Jul 2, 2021 at 14:23
  • Thank you for your attention and kind response, but when I tried this I got Error: [varnames<-,SpatRaster] cannot set these names
    – BeHana
    Commented Jul 2, 2021 at 14:41
  • After downloading your data, I too seem to have the same problem. My be best for now to store the names as a 'variable', like: rast.names <- names(test3), then assign these names to the raster names(LakeA) <- rast.names. Otherwise, seems to me an issue for the author to solve.
    – Stackbeans
    Commented Jul 2, 2021 at 14:52

1 Answer 1

2

You have a NetCDF file with many (52) variables (sub-datasets). When you open the file with rast these become "layers". Alternatively you can open the file with sds to keep the sub-dataset structure but that does not help you here (and you would need to skip the first two, see below).

library(terra)
f <- "ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190101-fv1.0.nc"
r <- rast(f)
r
#class       : SpatRaster 
#dimensions  : 21600, 43200, 52  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#sources     : ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190101-fv1.0.nc:water_surface_height_above_reference_datum  
              ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190101-fv1.0.nc:water_surface_height_uncertainty  
              ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190101-fv1.0.nc:lake_surface_water_extent  
              ... and 49 more source(s)
#varnames    : water_surface_height_above_reference_datum (water surface height above geoid) 
              water_surface_height_uncertainty (water surface height uncertainty) 
              lake_surface_water_extent (Lake Water Extent) 
              ...
#names       : water~datum, water~ainty, lake_~xtent, lake_~ainty, lake_~ature, lswt_~ainty, ... 
#unit        :           m,           m,         km2,         km2,      Kelvin,      Kelvin, ... 
#time        : 2019-01-01 

Note that there are 52 layers and sources (sub-datasets). There are names

head(names(r))
#[1] "water_surface_height_above_reference_datum" "water_surface_height_uncertainty"          
#[3] "lake_surface_water_extent"                  "lake_surface_water_extent_uncertainty"     
#[5] "lake_surface_water_temperature"             "lswt_uncertainty"                          

And also "longnames" (they are often much longer than the variable names, not in this case)

head(longnames(r))
# [1] "water surface height above geoid" "water surface height uncertainty" "Lake Water Extent"               
# [4] "Water extent uncertainty"         "lake surface skin temperature"    "Total uncertainty"               

You can also open the file with sds, but you need to skip "lon_bounds" and "lat_bounds" variables (dimensions)

s <- sds(f, 3:52)

Now read a vector data set (shapefile in this case) and crop

lake <- vect("hydro_p_LakeErie.shp")
rc <- crop(r, lake)
rc 

#class       : SpatRaster 
#dimensions  : 182, 555, 52  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : -83.475, -78.85, 41.38333, 42.9  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#names       : water~datum, water~ainty, lake_~xtent, lake_~ainty, lake_~ature, lswt_~ainty, ... 
#min values  :         NaN,         NaN,         NaN,         NaN,     271.170,       0.283, ... 
#max values  :         NaN,         NaN,         NaN,         NaN,     277.090,       0.622, ... 
#time        : 2019-01-01 
 

It can be convenient to save this to a GTiff file like this (or even better to use the filename argument in crop)

gtf <- writeRaster(rc, "test.tif", overwrite=TRUE)
gtf
#class       : SpatRaster 
#dimensions  : 182, 555, 52  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : -83.475, -78.85, 41.38333, 42.9  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : test.tif 
#names       : water~datum, water~ainty, lake_~xtent, lake_~ainty, lake_~ature, lswt_~ainty, ... 
#min values  :         NaN,         NaN,         NaN,         NaN,     271.170,       0.283, ... 
#max values  :         NaN,         NaN,         NaN,         NaN,     277.090,       0.622, ... 

What has changed is that the data are now in a file, rather then in memory. And you still have the layer (variable) names.

To write the layers as variables to a NetCDF file you need to create a SpatRasterDataset. You can do that like this:

x <- as.list(rc)
s <- sds(x)
names(s) <- names(rc)
longnames(s) <- longnames(r)
units(s) <- units(r)

Note the use of longnames(r) and units(r) (not rc). This is because r has subdatasets (and each has a longname and a unit) while rc does not.

Now use writeCDF

z <- writeCDF(s, "test.nc", overwrite=TRUE)
 
rc2 <- rast("test.nc")
rc2

#class       : SpatRaster 
#dimensions  : 182, 555, 52  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : -83.475, -78.85, 41.38333, 42.9  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#sources     : test.nc:water_surface_height_above_reference_datum  
              test.nc:water_surface_height_uncertainty  
              test.nc:lake_surface_water_extent  
              ... and 49 more source(s)
#varnames    : water_surface_height_above_reference_datum (water surface height above geoid) 
              water_surface_height_uncertainty (water surface height uncertainty) 
              lake_surface_water_extent (Lake Water Extent) 
              ...
#names       : water~datum, water~ainty, lake_~xtent, lake_~ainty, lake_~ature, lswt_~ainty, ... 
#unit        :           m,           m,         km2,         km2,      Kelvin,      Kelvin, ... 
#time        : 2019-01-01 

So it looks like we have a NetCDF with the same structure.

Note that the current CRAN version of terra drops the time variable if there is only one time step. The development version (1.3-11) keeps the time dimension, even of there is only one step.

You can install the development version with install.packages('terra', repos='https://rspatial.r-universe.dev')

3
  • Thank you. This helps to a certain level, but the time variable is important for me as I will be doing the clipping with several hundred files at different times so I need to retain the time variables.
    – BeHana
    Commented Jul 3, 2021 at 12:00
  • You can strore the date in the filename or use the development version (see updated answer) Commented Jul 3, 2021 at 18:12
  • Thank you very much! The development version worked! Now I will extend this cropping to all my files!
    – BeHana
    Commented Jul 4, 2021 at 10:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.