Flipped SpatRaster map when reading NetCDF file created by raster

[elgabbas]

Hello,
I want to read a netCDF file created by raster::writeRaster using terra::rast().

Create and save netCDF file

require(terra)
require(raster)

# creating raster
Map_raster <- raster::raster(
  nrows = 404, ncols = 390,  xmn = 2630000, xmx = 6530000,
  ymn = 1380000, ymx = 5420000, crs = "epsg:3035", vals = 1:(404*390)/1000)

# save as netcdf file
raster::writeRaster(
  Map_raster, "Map_raster.nc", overwrite=TRUE, format="CDF",  varname = "Map_raster")

Reading netCDF using raster

# as expected
plot(stack(Map_raster, raster("Map_raster.nc")), axes = F, box = F)

Reading netCDF using terra

Rast1 <- rast("Map_raster.nc")                  # Reading the NetCDF file using terra
crs(Rast1) <- "epsg:3035"
Rast2 <- rast(Map_raster)                         # coerce raster to SpatRaster 
plot(c(Rast1, Rast2))                                 # Rast1 is flipped!

# rast --> raster
# reading using terra::rast() then coerce to raster
plot(raster(rast("Map_raster.nc")))    # as expected!

# rast --> raster --> rast
plot(rast(raster(rast("Map_raster.nc"))))    # flipped again!

I am unsure if this issue is linked to this issue.

Your feedback is appreciated.

[elgabbas]

I can not really tell if the problem is with terra or raster.
I opened the resulting NetCDF file in QGIS and the European map is flipped!

[mablmabl]

A similar case with map flipped in y direction just using terra

download any file from this page for use in this example, I used April
https://data.marine.copernicus.eu/product/BALTICSEA_MULTIYEAR_PHY_003_011/files?path=BALTICSEA_MULTIYEAR_PHY_003_011%2Fcmems_mod_bal_phy_my_P1M-m_202303%2F2017%2F&subdataset=cmems_mod_bal_phy_my_P1D-m_202303

library(terra)  # I used version 1.7.74
library(ncdf4)  # I used version 1.22

file<-"BAL-MYP-NEMO_PHY-MonthlyMeans-201704.nc"

plot(terra::rast(file,subds = "sob"))               #gives correct map

nc_ds<-nc_open(file)
dim_lon <- ncvar_get(nc_ds, "lon")
dim_lat <- ncvar_get(nc_ds, "lat")
sob     <- ncvar_get(nc_ds, "sob")                  # salinity at bottom

plot(terra::rast(nrows=length(dim_lat), ncols=length(dim_lon), 
          xmin=min(dim_lon), xmax=max(dim_lon),
          ymin=min(dim_lat), ymax=max(dim_lat), 
          vals = as.vector(sob), 
          names = "sob",
          crs="EPSG:4326"))                         # gives map flipped in Y direction

Any clue why this happens?

[Rapsodia86]

Maybe this will help:
#1289

[mablmabl]

Many thanks for pointing me to #1289. Reading the nc file worked correctly with rast without any workaround. When I create a spatraster with terra the only workaround that worked was the last: flip vertical.

[rhijmans]

For what it is worth, @elgabbas' example seems to work with the lon/lat CRS.

library(terra)
library(raster)
r <- raster::raster(nrows = 404, ncols = 390, vals = 1:(404*390)/1000)
x <- raster::writeRaster(r, "map_raster.nc", overwrite=TRUE, format="CDF", varname = "Map_raster")
z <- rast("map_raster.nc")

It is unfortunate that it does not work correctly with other CRSs, but that would have to be fixed in "raster" (but that is not likely to happen anymore).

I do not really understand @mablmabl comment. That is, since

plot(terra::rast(file,subds = "sob")) 

"gives correct map"; terra would appear to be working as expected.