Netcdf files with unexpected order of time, easting and northing coordinates.

[ecor]

Hi all,
I would like to open a nc file and process it with terra package (see header below) . The terra::rast function (default settings) does not recognize the easting and northing coordinates but it takes "time coordinate as one of the the two. I see generally that the coordinate order is expected to be "time, easting, northing" and not "easting ,northing ,time" . Are there some specific settings in order terra::rast can read correctly the spatiotemporal gridded variables? Thank you for your attention and suggestion.

Best
Emanuele @ecor

$ ncdump -h METEO_DATA_202003120000.NC 
netcdf METEO_DATA_202003120000 {
dimensions:
        time = 24 ;
        northing = 180 ;
        easting = 185 ;
variables:
        double time(time) ;
                time:standard_name = "time" ;
                time:units = "days since 2020-03-12T00:00:00+01:00" ;
        double northing(northing) ;
                northing:units = "m" ;
                northing:long_name = "Northing" ;
                northing:coordinate\ system = "CH1903 / LV03 (EPSG:21781)" ;
        double easting(easting) ;
                easting:units = "m" ;
                easting:long_name = "Easting" ;
                easting:coordinate\ system = "CH1903 / LV03 (EPSG:21781)" ;
        double lwrc(easting, northing, time) ;
                lwrc:_FillValue = -9999. ;
                lwrc:standard_name = "downwelling_longwave_flux_in_air" ;
                lwrc:long_name = "Incoming longwave radiation" ;
                lwrc:coverage_content_type = "modelResult" ;
                lwrc:units = "W m-2" ;
        double sdri(easting, northing, time) ;
                sdri:_FillValue = -9999. ;
                sdri:standard_name = "direct_downwelling_shortwave_flux_in_air" ;
                sdri:long_name = "Incoming direct shortwave radiation on inclined surface" ;
                sdri:coverage_content_type = "modelResult" ;
                sdri:units = "W m-2" ;
        double sdfd(easting, northing, time) ;
                sdfd:_FillValue = -9999. ;
                sdfd:standard_name = "diffuse_downwelling_shortwave_flux_in_air" ;
                sdfd:long_name = "Incoming diffuse shortwave radiation" ;
                sdfd:coverage_content_type = "modelResult" ;
                sdfd:units = "W m-2" ;
        double tais(easting, northing, time) ;
                tais:_FillValue = -9999. ;
                tais:units = "K" ;
                tais:standard_name = "air_temperature" ;
                tais:long_name = "Air temperature at 10 m above ground" ;
                tais:coverage_content_type = "modelResult" ;
        double rhus(easting, northing, time) ;
                rhus:_FillValue = -9999. ;
                rhus:units = "%" ;
                rhus:standard_name = "relative_humidity" ;
                rhus:long_name = "Relative humidity at 10 m above ground" ;
                rhus:coverage_content_type = "modelResult" ;
        double pail(easting, northing, time) ;
                pail:_FillValue = -9999. ;
                pail:units = "Pa" ;
                pail:standard_name = "air_pressure" ;
                pail:long_name = "Air pressure" ;
                pail:coverage_content_type = "modelResult" ;
        double wnsd(easting, northing, time) ;
                wnsd:_FillValue = -9999. ;
                wnsd:units = "m/s" ;
                wnsd:standard_name = "wind_speed" ;
                wnsd:long_name = "Wind speed at 10 m above ground" ;
                wnsd:coverage_content_type = "modelResult" ;
        double wndd(easting, northing, time) ;
                wndd:_FillValue = -9999. ;
                wndd:units = "deg" ;
                wndd:standard_name = "wind_from_direction" ;
                wndd:long_name = "Wind direction at 10 m above ground" ;
                wndd:coverage_content_type = "modelResult" ;
        double proi(easting, northing, time) ;
                proi:_FillValue = -9999. ;
                proi:units = "mm/hr" ;
                proi:standard_name = "precipitation_flux" ;
                proi:long_name = "Precipitation" ;
                proi:coverage_content_type = "modelResult" ;

// global attributes:
                :title = "High-resolution hydrometeorological data for the Dischma valley and surrounding areas in eastern Switzerland" ;
                :summary = "This dataset contains high-resolution meteorological input data useful for developing and testing physics-based snow and landsurface models. The dataset contains hourly data and has a spatial resolution of 100m and an extent of 18500 by 18000m." ;
                :keywords = "downwelling_longwave_flux_in_air; direct_downwelling_shortwave_flux_in_air; diffuse_downwelling_shortwave_flux_in_air; air_temperature; relative_humidity; air_pressure; wind_speed; wind_from_direction; precipitation_flux" ;
                :conventions = "ACDD 1-3; CF-1.11" ;
                :creator_name = "Jan Magnusson; Tobias Jonas" ;
                :creator_email = "jan.magnusson@slf.ch; jonas@slf.ch" ;
                :creator_url = "https://orcid.org/0000-0003-0257-1862; https://orcid.org/0000-0003-0386-8676" ;
                :institution = "WSL Institute for Snow and Avalanche Research SLF" ;
                :project = "International Network for Alpine Research Catchment Hydrology" ;
                :geospatial_lat_min = 46.6734987688503 ;
                :geospatial_lat_max = 46.8353396574981 ;
                :geospatial_lon_min = 9.74259474964444 ;
                :geospatial_lon_max = 9.98461064390458 ;
                :date_created = "2024-08-08T23:22:05" ;
                :history = "File created at 2024-08-08T23:22:05 using Matlab by Jan Magnusson" ;
}
(r-reticulate)

[Breeze-Hu]

I may have had a similar problem at one point. rast doesn't recognise the time dimension very well, but rather is treated as longitude or latitude. At the time, my question on stack overflow was deleted by the forum administrator before it could be answered ....
One way I can think of is to extract the values, rotate the dimension, and then write the attributes in? (although I didn't try that subsequently)

[ecor]

Thanks, my question is about if rast is flexible and can open netcdf file in a non-standard format . In this case it takes "time"coordinate as latitude beacause "time" dimension is in the position expected for "latitude" dimension and so on. Otherwise variables can be directly extracted with "ncdf4" as array. Thanks @ecor

[Breeze-Hu]

I'm not sure I understood your point and question correctly, but I'll leave what I was thinking at the time. A 3D array of m*n*t (time) , I had problems recognising it as (m*t(time)*n),at the time I sought a way to expect rast to be able to manually define a certain dimension to be ‘time’ itself, but the question was removed.

Although it is offensive to mention other packs under terra, please forgive me. I ended up using xarray (a python package) for my processing.

[pvanlaake]

The ncdump listing clearly shows the nine data variables in the file and the three axes that they all use. So far, so good. The three axes (a.k.a. coordinate variable, composed of a dimension and a variable) are incompletely specified, however. The CF Metadata Conventions, which this file claims to adhere to, gives three options to identify how the coordinate variables are oriented in space, for X-Y: (1) the use of default units of "degrees_east" for longitude and "degrees_north" for latitude, which does not apply here because the data is in a projected CRS; (2) the use of a "standard_name" attribute, in this case with values "projected_x_coordinate" and "projected_y_coordinate"; or (3) the "axis" attribute with values "X" and "Y". The file provides neither of these mechanisms. The use of the "coordinate system" attribute is not standard and thus not recognised by a netCDF reader that applies the CF Metadata Conventions (and the name of the attribute is illegal too, having a space in it).

GDAL, which the terra package uses for reading netCDF files, understandably has trouble reading this file. It makes a number of guesses which are, unfortunately, not yielding the desired result. I downloaded the data you used (a link would have been useful, as well as notice on the web site that this was for a mere 31GB of data...) and this is what I get with terra::rast():

class       : SpatRaster 
dimensions  : 180, 24, 1665  (nrow, ncol, nlyr)
resolution  : 0.04166667, 100  (x, y)
extent      : -0.02083333, 0.9791667, 172000, 190000  (xmin, xmax, ymin, ymax)
coord. ref. :  
sources     : METEO_DATA_202003120000.NC:lwrc  (185 layers) 
              METEO_DATA_202003120000.NC:sdri  (185 layers) 
              METEO_DATA_202003120000.NC:sdfd  (185 layers) 
              ... and 6 more sources
varnames    : lwrc (Incoming longwave radiation) 
              sdri (Incoming direct shortwave radiation on inclined surface) 
              sdfd (Incoming diffuse shortwave radiation) 
              ...
names       : lwrc_~76050, lwrc_~76150, lwrc_~76250, lwrc_~76350, lwrc_~76450, lwrc_~76550, ... 
unit        :       W m-2,       W m-2,       W m-2,       W m-2,       W m-2,       W m-2, ... 
depth       : 776050 to 794450 (m [m]: 185 steps) 

The number of rows (northing) is correct but the columns are in fact the "time" axis and the number of layers is a merger of the easting axis and the 9 data variables in the file - the layer names show the result. The extent is also wrong and there is no CRS.

You can read the file more "cleanly" with package ncdfCF:

library(ncdfCF)

(ds <- open_ncdf("~/METEO_DATA_202003120000.NC"))
#> <Dataset> METEO_DATA_202003120000.NC 
#> Resource   : ~/METEO_DATA_202003120000.NC 
#> Format     : classic4 
#> Collection : Generic netCDF data 
#> Conventions: (not indicated) 
#> Keep open  : FALSE 
#> 
#> Variables:
#>  name long_name                                          units data_type  axes
#>  lwrc Incoming longwave radiation                        W m-2 NC_DOUBLE  time, northing, easting
#>  sdri Incoming direct shortwave radiation on inclined... W m-2 NC_DOUBLE  time, northing, easting
#>  sdfd Incoming diffuse shortwave radiation               W m-2 NC_DOUBLE  time, northing, easting
#>  tais Air temperature at 10 m above ground               K     NC_DOUBLE  time, northing, easting
#>  rhus Relative humidity at 10 m above ground             %     NC_DOUBLE  time, northing, easting
#>  pail Air pressure                                       Pa    NC_DOUBLE  time, northing, easting
#>  wnsd Wind speed at 10 m above ground                    m/s   NC_DOUBLE  time, northing, easting
#>  wndd Wind direction at 10 m above ground                deg   NC_DOUBLE  time, northing, easting
#>  proi Precipitation                                      mm/hr NC_DOUBLE  time, northing, easting
#> 
#> Attributes:
#>  name               type      length value 
#>  title              NC_CHAR   108    High-resolution hydrometeorological data for th...
#>  summary            NC_CHAR   244    This dataset contains high-resolution meteorolo...
#>  keywords           NC_CHAR   220    downwelling_longwave_flux_in_air; direct_downwe...
#>  conventions        NC_CHAR    17    ACDD 1-3; CF-1.11
#>  creator_name       NC_CHAR    27    Jan Magnusson; Tobias Jonas
#>  creator_email      NC_CHAR    34    jan.magnusson@slf.ch; jonas@slf.ch
#>  creator_url        NC_CHAR    76    https://orcid.org/0000-0003-0257-1862; https://...
#>  institution        NC_CHAR    49    WSL Institute for Snow and Avalanche Research SLF 
#>  project            NC_CHAR    61    International Network for Alpine Research Catch...
#>  geospatial_lat_min NC_DOUBLE   1    46.6734987688503 
#>  geospatial_lat_max NC_DOUBLE   1    46.8353396574981
#>  geospatial_lon_min NC_DOUBLE   1    9.74259474964444
#>  geospatial_lon_max NC_DOUBLE   1    9.98461064390458 
#>  date_created       NC_CHAR    19    2024-08-08T23:22:05
#>  history            NC_CHAR    65    File created at 2024-08-08T23:22:05 using Matla...

# Open a data variable
(lwrc <- ds[["lwrc"]])
#> <Variable> lwrc 
#> Long name: Incoming longwave radiation 
#> 
#> Axes:
#>  axis name     long_name length values                                        unit
#>  T    time                24    [2020-03-12T00:00:00 ... 2020-03-12T23:00:00] days since 2020-03-12T00:00:00+01:00
#>       northing Northing  180    [189950 ... 172050]                           m
#>       easting  Easting   185    [776050 ... 794450]                           m
#> 
#> Attributes:
#>  name                  type      length value                           
#>  _FillValue            NC_DOUBLE  1     -9999                           
#>  standard_name         NC_CHAR   32     downwelling_longwave_flux_in_air
#>  long_name             NC_CHAR   27     Incoming longwave radiation     
#>  coverage_content_type NC_CHAR   11     modelResult                     
#>  units                 NC_CHAR    5     W m-2

This all looks ok but the northing and easting axes are not recognised as Y and X. You need to remedy that manually:

lwrc$axes[["northing"]]$orientation <- "Y"
lwrc$axes[["easting"]]$orientation <- "X"
lwrc
#> <Variable> lwrc 
#> Long name: Incoming longwave radiation 
#> 
#> Axes:
#>  axis name     long_name length values                                        unit
#>  T    time                24    [2020-03-12T00:00:00 ... 2020-03-12T23:00:00] days since 2020-03-12T00:00:00+01:00
#>  Y    northing Northing  180    [189950 ... 172050]                           m
#>  X    easting  Easting   185    [776050 ... 794450]                           m
#> 
#> Attributes:
#>  name                  type      length value                           
#>  _FillValue            NC_DOUBLE  1     -9999                           
#>  standard_name         NC_CHAR   32     downwelling_longwave_flux_in_air
#>  long_name             NC_CHAR   27     Incoming longwave radiation     
#>  coverage_content_type NC_CHAR   11     modelResult                     
#>  units                 NC_CHAR    5     W m-2

(Note that you have to do this only once: all nine data variables share the same instance of northing and easting and it is those objects that the $orentation applies to.)

Now you have a fully specified coordinate system. You can easily convert this data variable to a terra::SpatRaster:

(r <- lwrc$data()$terra())
#> class       : SpatRaster 
#> dimensions  : 180, 185, 24  (nrow, ncol, nlyr)
#> resolution  : 100, 100  (x, y)
#> extent      : 776000, 794500, 172000, 190000  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> names       : 2020-~00:00, 2020-~00:00, 2020-~00:00, 2020-~00:00, 2020-~00:00, 2020-~00:00, ... 
#> min values  :       173.8,       174.1,       174.4,       173.1,       177.0,       180.6, ... 
#> max values  :       223.3,       222.6,       221.5,       220.1,       221.8,       223.7, ...

All done, except for the CRS. Given the extent of the data variable, it is a bit odd that EPSG:4326 is assigned. You can set your own CRS with r$crs <- "epsg:21781" but it appears that terra does not recognise that even though it is a valid EPSG code.

In conclusion, this is an issue with the data file and not with terra.

[mdsumner]

in an upcoming GDAL release (3.12, a while until that's available generally) you might get away with

vrt://METEO_DATA_202003120000.NC?transpose=lwrc:2,1

as the input to rast() rather than the file name itself. With earlier GDAL you can use also use cli to do this with

gdalmdimtranslate in.nc out.nc -array "name=lwrc,transpose=[2,1,0]"

I suspect it would work better now if the file had a "grid_mapping" var for the crs. Can you point to this file so I can check this out? (I wonder if GDAL could benefit from a virtualization of 'grid_mapping', so we could inject what it should use via VRT)

(I found, I think, files here: https://www.envidat.ch/#/metadata/inarch-dischma-dataset)

[ecor]

Thanks for the feedback. Yes they comes from Dischma catchment meteorological dataset. At the end working directly within an R session, I opened each using directly ncdf4 directly, modified the attributes and the variables (see lines of code below) and then created a new netCDF file for each original one. The new one can be obened by _terra::rast_ as a raster map. I share the solution here:

#
# DISCHMA METEO DATASET 
#
library(ncdf4)
library(terra)
library(stringr)
library(lubridate)
### DISCHMA DATASET: https://www.envidat.ch/#/metadata/inarch-dischma-dataset


meteo_dir <- ".../dischma/meteo/"
output_file_path <- ".../dischma_new/meteo/%s" ## C formatted string with %s

meteo_files <- meteo_dir |> list.files(pattern=".NC",full.names=TRUE)

## try the fist ten nc files: 
meteo_files <- meteo_files[1:10] ## you can comment for all files 

for (meteo_file in meteo_files) {
# Input and output file names
input_file <- meteo_file ##"METEO_DATA_202003120000.NC"
output_file <- input_file |> basename() |> sprintf(fmt=output_file_path) ##"METEO_DATA_202003120000_reordered.nc"

# Open the original NetCDF file
nc_in <- nc_open(input_file)

# Read dimension values
time_vals <- ncvar_get(nc_in, "time")
northing_vals <- ncvar_get(nc_in, "northing")
easting_vals <- ncvar_get(nc_in, "easting")

# Define new dimensions in the desired order
dim_time <- ncdim_def("time", ncatt_get(nc_in, "time", "units")$value, time_vals)
dim_easting <- ncdim_def("easting", ncatt_get(nc_in, "easting", "units")$value, easting_vals)
dim_northing <- ncdim_def("northing", ncatt_get(nc_in, "northing", "units")$value, northing_vals)

# List of variables to reorder
vars_to_reorder <- c("lwrc", "sdri", "sdfd", "tais", "rhus", "pail", "wnsd", "wndd", "proi")

# Define variables with reordered dimensions
var_defs <- lapply(vars_to_reorder, function(varname) {
  var <- nc_in$var[[varname]]
  ncvar_def(name = var$name,
            units = ncatt_get(nc_in, varname, "units")$value,
            dim = list(dim_easting, dim_northing,dim_time), ### EC 20250519
            missval = var$missval,
            prec = "double")
})

var_defs$crs  <- ncvar_def(
  name = "crs",
  units = "",
  dim = list(),  # NUll dimension = scalare variable
  longname = "Coordinate Reference System",
  prec = "integer"
)

#ncvar_put(nc_out,crs_var, 1)
## Create new NetCDF file
nc_out <- nc_create(output_file, var_defs,force_v4=FALSE)

# Write dimension variables
ncvar_put(nc_out, "time", time_vals)
ncvar_put(nc_out, "easting", easting_vals)
ncvar_put(nc_out, "northing", northing_vals)

# Copy dimension attributes
for (dim_name in c("time", "easting", "northing")) {
  att_names <- names(ncatt_get(nc_in, dim_name))
  for (att in att_names) {
    ncatt_put(nc_out, dim_name, att, ncatt_get(nc_in, dim_name, att)$value)
    ## PUT CRS
   
    
  }
  if (dim_name %in% c("northing","easting")) {
     
     attrcrs <- ncatt_get(nc_in, dim_name, att)$value
    
     epsg_code <- sub(".*\\(EPSG:(\\d+)\\).*", "\\1", attrcrs) |> as.integer()
     ###ncatt_put(nc_out,"crs","code" , epsg_code)
     wkt <- sf::st_crs(epsg_code)$wkt
     ncatt_put(nc_out,"crs","crs_wkt",wkt)
   } else if (dim_name=="time") {
     attrunits_time <- ncatt_get(nc_in, dim_name,"units")$value
     t0 <- attrunits_time |> str_split("since") |> sapply(FUN=function(x){x[2]}) |> 
       str_trim() |> str_replace_all(":","") |> 
       as.POSIXct(format="%Y-%m-%dT%H%M%S%z",tz="UTC")
     time_vals2 <- t0+seconds(round(time_vals*24*60*60))
     t00 <- as.Date(time_vals2[3]) |> as.POSIXct(tz="UTC")
     time_vals3 <- difftime(time_vals2,t00,units="secs") |> as.numeric()
     ncvar_put(nc_out, "time", time_vals3)
     
     
     attrunits_time2 <- "seconds since %s" |> sprintf(as.character(as.Date(t00)))
     ncatt_put(nc_out, dim_name,"units",attrunits_time2)
   }   
  
}

ncatt_put(nc_out, "easting", "standard_name", "projection_x_coordinate")
ncatt_put(nc_out, "northing", "standard_name", "projection_y_coordinate")

# Write variables and copy all attributes
for (varname in vars_to_reorder) {
  data <- ncvar_get(nc_in, varname)
  data_reordered <- aperm(data, c(3, 2, 1))  # from (easting 3 , northing 2 , time 1 ) to (time 1 , easting 3 , northing2 )
  ncvar_put(nc_out, varname, data_reordered)
  
  # Copy all variable attributes
  att_names <- names(ncatt_get(nc_in, varname))
  for (att in att_names) {
    ncatt_put(nc_out, varname, att, ncatt_get(nc_in, varname, att)$value)
    
  }
 
  ncatt_put(nc_out, varname,"grid_mapping","crs")
}

# Copy global attributes
global_atts <- names(ncatt_get(nc_in, 0))
for (att in global_atts) {
  ncatt_put(nc_out, 0, att, ncatt_get(nc_in, 0, att)$value)
}


# Close files
nc_close(nc_in)
nc_close(nc_out)

cat("✅ NetCDF file successfully reordered to (time, easting, northing) and saved as:", output_file, "\n")
}

and then it works:

> output_file
[1] ".../dischma_new/meteo/METEO_DATA_201610100000.NC"
> rast(output_file)
class       : SpatRaster 
size        : 180, 185, 216  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 776000, 794500, 172000, 190000  (xmin, xmax, ymin, ymax)
coord. ref. : CH1903 / LV03 (EPSG:21781) 
sources     : METEO_DATA_201610100000.NC:lwrc  (24 layers) 
              METEO_DATA_201610100000.NC:sdri  (24 layers) 
              METEO_DATA_201610100000.NC:sdfd  (24 layers) 
              ... and 6 more sources
varnames    : lwrc (Incoming longwave radiation) 
              sdri (Incoming direct shortwave radiation on inclined surface) 
              sdfd (Incoming diffuse shortwave radiation) 
              ...
names       : lwrc_1, lwrc_2, lwrc_3, lwrc_4, lwrc_5, lwrc_6, ... 
unit        :  W m-2,  W m-2,  W m-2,  W m-2,  W m-2,  W m-2, ... 
time        : 2016-10-09 23:00:00 to 2016-10-10 22:00:00 UTC (24 steps) 

[ecor]

This is not a direct and internal issue of terra but I hope it can be useful in case on how to use some non-standard netCDF files. Thank you
best @ecor