writeCDF SpatRaster names

[ForChimneySwifts]

When writing a SpatRasterDataset to a .nc file, could the writeCDF function include the original layer names within the new layer names created by the function?
For example, currently, layer names are renamed by the writeCDF function to be varname_1, varname_2, etc. Instead of using this convention, could this function rename these layers to varname_"original layername" instead, where "original_layername" is the layer's name in the SDS object? It is frustrating to have to break up an SDS into many spatRasters just so that I can save the layer names correctly in the nc file.
Or is there a way to do this with the writeCDF function that I am unaware of?
Thank you

[rhijmans]

Can you include a small example that illustrates what you are seeing?

[ForChimneySwifts]

I can upload example files, but I thought this screenshot might be enough. As you can see, the names of the TH layers (elevations) have been changed from "Elev=0", "Elev=0.2", etc to generic names after saving as a netcdf file.

I know I can use the "time" variable to save unique datetime attributes for each layer, but there is no way that I can see to save unique elevation information for each layer. This information is lost when converting to a nc using writeCDF.

One solution would be if writeCDF could append the original name attribute to the generic "name" being created.

example.pdf

[rhijmans]

I would like a concrete example, either generated with code or with a file. We need to be sure that we talk about the same things. ncdf files have variable names, not layer names. Perhaps the issue is that terra needs to understand higher dimensions, but I am not sure what it is without a clear example, so that we have a well defined behavior, and clear expectations.

[rhijmans]

Also, please never include screenshots in an issue (except for a plot, where relevant). You can copy and paste the exact input/output from the console instead.

[ForChimneySwifts]

library(terra)
m <- matrix(1:25, nrow=5, ncol=5)
rm <- rast(m)
rm <- c(rm,rm,rm)
names(rm) <- c("elev_0","elev_0.5","elev_0.7")
test_rm  <- sds(rm,rm,rm,rm)
names(test_rm) <- c("TH","DBZH","VRADH","filtered")
writeCDF(test_rm, filename = "test_rm.nc",overwrite=T)
test_nc <- sds("test_rm.nc")

test_rm$TH
#class : SpatRaster
#dimensions : 5, 5, 3 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : 0, 5, 0, 5 (xmin, xmax, ymin, ymax)
#coord. ref. :
#source(s) : memory
#names : elev_0, elev_0.5, elev_0.7
#min values : 1, 1, 1
#max values : 25, 25, 25

test_nc$TH
#class : SpatRaster
#dimensions : 5, 5, 3 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : 0, 5, 0, 5 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
#source : test_rm.nc:TH
#varname : TH
#names : TH_1, TH_2, TH_3

[rhijmans]

In this example you have variables c("TH","DBZH","VRADH","filtered") . You can change the first variable to "elev" to get

names(test_nc[1])[1:3]
# names       : elev_1, elev_2, elev_3 

The layer names are derived from the variable name of the array. The NCDF CF-1standard does not support "layer names" (names for the third dimension of an x, y, z array). It would be possible to make {terra} write these and read them again.

Perhaps a better approach is to automatically split SpatRasters based on their layer names. Either in writeCDF or with a function you can call before that. That could keep pattern_1 .. pattern_n together as variable pattern. In your case you would get 12 subdatsets: TH_elev_0, TH_elev_0.5, TH_elev_0.7, DBZH_elev_0, DBZH_elev_0.5, DBZH_elev_0.7, VRADH_elev_0, VRADH_elev_0.5, VRADH_elev_0.7, filtered_elev_0, filtered_elev_0.5, filtered_elev_0.7.

And you can get that with current terra by writing a SpatRaster instead of a SpatRasterDataset and using split=TRUE :

writeCDF(rast(test_rm), filename = "test2_rm.nc",overwrite=T, split=TRUE)
rast("test2_rm.nc")
#class       : SpatRaster 
#dimensions  : 5, 5, 12  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 5, 0, 5  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#sources     : test2_rm.nc:TH_elev_0  
#              test2_rm.nc:TH_elev_0.5  
#              test2_rm.nc:TH_elev_0.7  
#              ... and 9 more sources
#varnames    : TH_elev_0 
#              TH_elev_0.5 
#              TH_elev_0.7 
#              ...
#names       : TH_elev_0, TH_elev_0.5, TH_elev_0.7, DBZH_elev_0, DBZH_~v_0.5, DBZH_~v_0.7, ... 

[ForChimneySwifts]

Is "time" not a third dimension of an x,y,z array? How is time stored, and could it not be generic?

[rhijmans]

time is a dimension, like x/longitude and y/latitude. There are no names associated with each step. There is not a name for each latitude step. The third dimension could also be, for example, depth and we would have the depth values (10, 20), but not names like "shallow" and "deep"

[rhijmans]

Another thing you can do is change the label "time" to another label to get the below:

(writeCDF(rm, "test.nc", zname="depth", overwrite=T))
#class       : SpatRaster 
#dimensions  : 5, 5, 3  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 5, 0, 5  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source      : test.nc 
#names       : test_depth=1, test_depth=2, test_depth=3 

[ForChimneySwifts]

Thank you for your explanations. I only have one further question ... Is it possible to use depths (or angles, in my example) as the third (time) dimension in terra? Are the unique values of the 3rd dimension stored in the netcdf, or only the start value and step value?

[rhijmans]

The unique values are stored. Right now only time is used (but that I can fix). You can use time as long as your values are integers

time(rm) <- c(10, 20, 30)
(writeCDF(rm, "test.nc", zname="angle", overwrite=T))
#class       : SpatRaster 
#dimensions  : 5, 5, 3  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 5, 0, 5  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source      : test.nc 
#names       : test_angle=10, test_angle=20, test_angle=30 

But it would be better to allow doing this via "depth" because it allows real numbers; allow setting the "depth-name" to e.g., "angle", and reading the file again. There could also be other options for creating the layer names.

In your example you would name the third dimension "elev" and you would choose not to prepend the variable name (that should not be necessary anyway if you read the file as an SpatRasterDataset.

And it should be possible to use depth and time to write four-dimensional data.

[ForChimneySwifts]

Would the 3rd dimension always have to be an integer? Or could it be any value, as long as these values were defined in the "time" variable? For example
elevation(rm) <- c(-0.5, -0.3, 0, 0.3, 0.5)

[rhijmans]

With the last commit, a small step forward, depth is used if there is no time dimension; allowing the use of real numbers:

library(terra)
rm <- rast(nrow=5, ncol=5, vals=1:75, nlyr=3)
depth(rm) <- c(0, 0.5, 0.7)
names(rm) <- c("elev_0","elev_0.5","elev_0.7")
srm  <- sds(rm,rm,rm,rm)
names(srm) <- c("TH","DBZH","VRADH","filtered")
writeCDF(srm, filename = "test.nc", zname="angle", overwrite=T)
s <- sds("test.nc")
s[1]

#class       : SpatRaster 
#dimensions  : 5, 5, 3  (nrow, ncol, nlyr)
#resolution  : 72, 36  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source      : test.nc:TH 
#varname     : TH 
#names       : TH_angle=0, TH_angle=0.5, TH_angle=0.7 

[ForChimneySwifts]

Great! Could each dataset within an SDS have a different third dimension?

depth(rm) <- c(0,0.5,0.7)
depth(rm1) <- c(-0.5, 0)
srm  <- sds(rm,rm1)
names(srm) <- c("TH","DBZH")
writeCDF(srm, filename = "test.nc", zname="angle", overwrite=T)
s <- sds("test.nc")

[rhijmans]

depth is the generic name for the third dimension, and you can now set a specific name for that dimension, that is used when writing a ncdf file. time is the fourth dimension (or third if there is no depth dimension) and its name can be changed with zname. Some more work is needed, but I think this is moving in the right direction.

library(terra)
r1 <- rast(nrow=5, ncol=5, vals=1:75, nlyr=3)
depth(r1) <- c(0, 0.5, 0.7)
r2 <- rast(nrow=5, ncol=5, vals=1:50, nlyr=2)
depth(r2) <- c(-10, 20)
depthName(r1) <- "angle"
depthName(r2) <- "height"
srm <- sds(r1, r2)
names(srm) <- c("TH", "DBZH")
x <- writeCDF(srm, filename = "test.nc", overwrite=T)

x[1]
#class       : SpatRaster 
#dimensions  : 5, 5, 3  (nrow, ncol, nlyr)
#resolution  : 72, 36  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source      : test.nc:TH 
#varname     : TH 
#names       : TH_angle=0, TH_angle=0.5, TH_angle=0.7 
#angle       : 0 to 0.7 

x[2]
#class       : SpatRaster 
#dimensions  : 5, 5, 2  (nrow, ncol, nlyr)
#resolution  : 72, 36  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source      : test.nc:DBZH 
#varname     : DBZH 
#names       : DBZH_height=-10, DBZH_height=20 
#height      : -10 to 20 

[ForChimneySwifts]

Thank you!

[rhijmans]

It would be helpful to know more about what you need. Perhaps terra can support that better. RH

…

On Mon, Mar 3, 2025 at 6:47 PM Breeze-Hu ***@***.***> wrote: I used to get confused on layer names for high dimensional data as well. Now one of my solution is to use python xarray for data reading open_dataset ,and then do extraction for dimension names. Then use terra::names function for corresponding naming. And in using terra I try to use 3D data in order to avoid my own dimensional problems. — Reply to this email directly, view it on GitHub <#1756 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACXTCNJSI3J5NMKFVPWWK632SUH55AVCNFSM6AAAAABYH5YDGGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMOJWGA2DANBSGM> . You are receiving this because you commented.Message ID: ***@***.***> [image: Breeze-Hu]*Breeze-Hu* left a comment (rspatial/terra#1756) <#1756 (comment)> I used to get confused on layer names for high dimensional data as well. Now one of my solution is to use python xarray for data reading open_dataset ,and then do extraction for dimension names. Then use terra::names function for corresponding naming. And in using terra I try to use 3D data in order to avoid my own dimensional problems. — Reply to this email directly, view it on GitHub <#1756 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACXTCNJSI3J5NMKFVPWWK632SUH55AVCNFSM6AAAAABYH5YDGGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMOJWGA2DANBSGM> . You are receiving this because you commented.Message ID: ***@***.***>