Extent and crs of rasters not written by writeCDF

user_4003696 0
Viewed 233

When writing rasters in netCDF files, I always get the warning message: "[rast] unknown extent". Indeed, the extent is not written in the external file. Neither is the crs.

library(terra)
#terra version 1.0.2

r <- rast(ncol=2, nrow=2, vals=c(5.3, 7.1, 3, 1.2))
crs(r)<-"epsg:27572"
ext(r)
#SpatExtent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

t<-writeCDF(r,"test.ncdf",overwrite=TRUE)
#Warning message:
#[rast] unknown extent
 
ext(t)  # extension is not correct
#SpatExtent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)

crs(t)  # crs is not correct
#[1] "GEOGCRS[\"unknown\",\n    DATUM[\"World Geodetic System 1984\",\n     ...

Perhaps there is a peculiar syntax to use here. I explored ?writeCDF, but could not find any clue.

terra raster netcdf
1 Answers

This points at an issue with GDAL --- depending on whether you think that .ncdf is a common filename extension for netCDF files.

library(terra)
#terra version 1.0.3
r <- rast(ncol=2, nrow=2, vals=c(5.3, 7.1, 3, 1.2))

Note the different file extensions, .nc, .cdf, .ncdf or missing.

# ok
x <- writeCDF(r, "test1.nc", overwrite=TRUE)
y <- writeCDF(r, "test2.cdf", overwrite=TRUE)

# not ok
z <- writeCDF(r, "test3.ncdf", overwrite=TRUE)
#Warning message:
#[rast] unknown extent
a <- writeCDF(r, "test4", overwrite=TRUE)
#Warning message:
#[rast] unknown extent

GDALinfo shows:

describe("test1.nc")[1] 
#[1] "Driver: netCDF/Network Common Data Format"
describe("test3.ncdf")[1]
#[1] "Driver: HDF5Image/HDF5 Dataset"

It looks like GDAL first tries the netCDF driver when the extension is .nc or .cdf, but that it first tries the HDF5 driver when it is .ncdf or missing --- and since this does not fail (the warning comes from terra, not from GDAL), that is what it uses.

This is the GDAL version on windows.

gdal()
#[1] "3.0.4"

I see the same behavior with GDAL 2.2.3 on linux and 3.2.0 on mac.

You can work around that by not using .ncdf or by specifying the driver when opening the file:

rast('NETCDF:"test3.ncdf"')
#class       : SpatRaster 
#dimensions  : 2, 2, 1  (nrow, ncol, nlyr)
#resolution  : 180, 90  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : NETCDF:test1.ncdf 
#varname     : test1 
#name        : test1

I do not think there is anything wrong with the CRS (it is the same as crs(r)). However, I should note that terra writes the proj4 and wkt strings to a ncdf file, and does not follow the ncdf standard in that respect.

(You are asking a question about a method that is only available in the development version of terra. I appreciated that very much, but raising an issue on the terra github site would be more appropriate in this case. I will make writeCDF give a warning when the file extension is not .nc or .cdf)

user_6352450
Related