This NetCDF has a spherical CRS longlat crs in "crs_wkt", but the NetCDF driver hardcodes WGS84 as the geolocation SRS:
gdalinfo /vsicurl/https://github.com/rgdal-dev/rasterwise/raw/master/data-raw/spherical.nc
...
Geolocation:
SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
X
...I believe it is hardcoded here: https://github.com/OSGeo/gdal/blob/35c434527aab4acb86af66480db00d4fec4b5400/frmts/netcdf/netcdfdataset.cpp#L4805
spherical.nc created with this:
(note the data in gebco var is not oriented correctly but that doesn't affect the purpose of this example, I will fix at some point)
dm <- c(128, 128)
d <- vapour::gdal_raster_data(sds::gebco(), target_dim = c(128L, 128L), target_crs = "+proj=laea")
#> Warning in warp_general_cpp(dsn, target_crs, as.numeric(target_ext),
#> as.integer(target_dim), : no source crs, target crs is ignored
crs <- vapour::vapour_srs_wkt("+proj=longlat +R=6378137")
library(RNetCDF)
nc <- create.nc("spherical.nc", format = "netcdf4")
x <- dim.def.nc(nc, "x", dm[1L])
y <- dim.def.nc(nc, "y", dm[2L])
gebco <- var.def.nc(nc, "gebco", "NC_FLOAT", c("x", "y"))
att.put.nc(nc, gebco, "grid_mapping", "NC_STRING", "crs")
att.put.nc(nc, gebco, "coordinates", "NC_STRING", "lon lat")
crsvar <- "crs"
var <- var.def.nc(nc, crsvar, "NC_INT", NA)
var.put.nc(nc, crsvar, 1L)
att.put.nc(nc, crsvar, "comment", "NC_STRING",
"This is a container variable that describes the grid_mapping used by the data in this file. This variable does not contain any data; only information about the geographic coordinate system."
)
att.put.nc(nc, crsvar, "grid_mapping_name", "NC_STRING", "latitude_longitude")
att.put.nc(nc, crsvar, "crs_wkt", "NC_STRING", crs)
xx <- var.def.nc(nc, "lon", "NC_DOUBLE", c("x", "y"))
yy <- var.def.nc(nc, "lat", "NC_DOUBLE", c("x", "y"))
att.put.nc(nc, xx, "long_name", "NC_STRING", "longitude")
att.put.nc(nc, yy, "long_name", "NC_STRING", "latitude")
xycoords <- vaster::xy(attr(d, "dimension"), attr(d, "extent"))
xycoords <- reproj::reproj_xy(xycoords, crs, source = attr(d, "projection"))
#> Error detected, some values Inf (error code: 2050)
#>
#> ' Point outside of projection domain
#>
#> '
var.put.nc(nc, xx, matrix(xycoords[,1L], dm[2L], byrow = TRUE))
var.put.nc(nc, yy, matrix(xycoords[,2L], dm[2L], byrow = TRUE))
var.put.nc(nc, gebco, matrix(d[[1L]], dm[2L], byrow = TRUE))
close.nc(nc)
writeLines(system("ncdump -h spherical.nc", intern = TRUE))
#> netcdf spherical {
#> dimensions:
#> x = 128 ;
#> y = 128 ;
#> variables:
#> float gebco(y, x) ;
#> string gebco:grid_mapping = "crs" ;
#> string gebco:coordinates = "lon lat" ;
#> int crs ;
#> string crs:comment = "This is a container variable that describes the grid_mapping used by the data in this file. This variable does not contain any data; only information about the geographic coordinate system." ;
#> string crs:grid_mapping_name = "latitude_longitude" ;
#> string crs:crs_wkt = "GEOGCS[\"unknown\",DATUM[\"unknown\",SPHEROID[\"unknown\",6378137,0]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Longitude\",EAST],AXIS[\"Latitude\",NORTH]]" ;
#> double lon(y, x) ;
#> string lon:long_name = "longitude" ;
#> double lat(y, x) ;
#> string lat:long_name = "latitude" ;
#> }
writeLines(system("gdalinfo NETCDF:spherical.nc:gebco", intern = TRUE))
#> Driver: netCDF/Network Common Data Format
#> Files: spherical.nc
#> Size is 128, 128
#> Coordinate System is:
#> GEOGCRS["unknown",
#> DATUM["unknown",
#> ELLIPSOID["unknown",6378137,0,
#> LENGTHUNIT["metre",1,
#> ID["EPSG",9001]]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]],
#> CS[ellipsoidal,2],
#> AXIS["longitude",east,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> AXIS["latitude",north,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]]]
#> Data axis to CRS axis mapping: 1,2
#> Metadata:
#> crs#comment=This is a container variable that describes the grid_mapping used by the data in this file. This variable does not contain any data; only information about the geographic coordinate system.
#> crs#crs_wkt=GEOGCS["unknown",DATUM["unknown",SPHEROID["unknown",6378137,0]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]
#> crs#grid_mapping_name=latitude_longitude
#> gebco#coordinates=lon lat
#> gebco#grid_mapping=crs
#> Geolocation:
#> LINE_OFFSET=0
#> LINE_STEP=1
#> PIXEL_OFFSET=0
#> PIXEL_STEP=1
#> SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
#> X_BAND=1
#> X_DATASET=NETCDF:"spherical.nc":lon
#> Y_BAND=1
#> Y_DATASET=NETCDF:"spherical.nc":lat
#> Corner Coordinates:
#> Upper Left ( 0.0, 0.0)
#> Lower Left ( 0.0, 128.0)
#> Upper Right ( 128.0, 0.0)
#> Lower Right ( 128.0, 128.0)
#> Center ( 64.0, 64.0)
#> Band 1 Block=128x1 Type=Float32, ColorInterp=Undefined
#> NoData Value=9.96920996838686905e+36
#> Metadata:
#> coordinates=lon lat
#> grid_mapping=crs
#> NETCDF_VARNAME=gebcoCreated on 2024-03-22 with reprex v2.0.2