I'm working with data that comes in a deprecated crs and want to avoid reprojection after converting the dataset to a SpatRaster. It seems that gdal is automatically replacing EPSG:2163 with EPSG:9311 (see here). Gdal apparently displays a warning, but the warning is not displayed by terra.
I can force the SpatRaster into the correct crs by using the full crs string. Subsequent operations (e.g., mask) appear to use the deprecated crs, as desired.
r <- rast(ncol = 1, nrow = 1, vals = 1,
xmin = -1694992.5, xmax = -1694137.5,
ymin = -430492.5, ymax = -429367.5,
crs = 'EPSG:2163') # creates raster in EPSG:9311
crs(r) <- 'EPSG:2163' # doesn't work
# however, this successfully forces crs to EPSG:2163:
crs(r) <-"PROJCRS[\"US National Atlas Equal Area\",\n BASEGEOGCRS[\"Unspecified datum based upon the Clarke 1866 Authalic Sphere\",\n DATUM[\"Not specified (based on Clarke 1866 Authalic Sphere)\",\n ELLIPSOID[\"Clarke 1866 Authalic Sphere\",6370997,0,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4052]],\n CONVERSION[\"US National Atlas Equal Area\",\n METHOD[\"Lambert Azimuthal Equal Area (Spherical)\",\n ID[\"EPSG\",1027]],\n PARAMETER[\"Latitude of natural origin\",45,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-100,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"False easting\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"easting (X)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (Y)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Statistical analysis.\"],\n AREA[\"United States (USA) - onshore and offshore.\"],\n BBOX[15.56,167.65,74.71,-65.69]],\n ID[\"EPSG\",2163]]"
However, when I write the raster to a GeoTIFF, I haven't yet figured out how to override gdal replacing the crs.
tf <- tempfile(fileext = '.tif')
writeRaster(r, tf) # this doesn't work
writeRaster(r, tf, gdal = 'OSR_USE_NON_DEPRECATED=NO') # neither does this
The clunky solution is to just use the full crs string to redefine the crs after loading the GeoTiff into R. Is there a better way to override crs replacement for both writing and reading the file?
You example data
library(terra)
r <- rast(ncol = 1, nrow = 1, vals = 1,
xmin = -1694992.5, xmax = -1694137.5,
ymin = -430492.5, ymax = -429367.5,
crs = 'EPSG:2163') # creates raster in EPSG:9311
crs(r) <-"PROJCRS[\"US National Atlas Equal Area\", BASEGEOGCRS[\"Unspecified datum based upon the Clarke 1866 Authalic Sphere\", DATUM[\"Not specified (based on Clarke 1866 Authalic Sphere)\", ELLIPSOID[\"Clarke 1866 Authalic Sphere\",6370997,0, LENGTHUNIT[\"metre\",1]]], PRIMEM[\"Greenwich\",0, ANGLEUNIT[\"degree\",0.0174532925199433]], ID[\"EPSG\",4052]], CONVERSION[\"US National Atlas Equal Area\", METHOD[\"Lambert Azimuthal Equal Area (Spherical)\", ID[\"EPSG\",1027]], PARAMETER[\"Latitude of natural origin\",45, ANGLEUNIT[\"degree\",0.0174532925199433], ID[\"EPSG\",8801]], PARAMETER[\"Longitude of natural origin\",-100, ANGLEUNIT[\"degree\",0.0174532925199433], ID[\"EPSG\",8802]], PARAMETER[\"False easting\",0, LENGTHUNIT[\"metre\",1], ID[\"EPSG\",8806]], PARAMETER[\"False northing\",0, LENGTHUNIT[\"metre\",1], ID[\"EPSG\",8807]]], CS[Cartesian,2], AXIS[\"easting (X)\",east, ORDER[1], LENGTHUNIT[\"metre\",1]], AXIS[\"northing (Y)\",north, ORDER[2], LENGTHUNIT[\"metre\",1]], USAGE[ SCOPE[\"Statistical analysis.\"], AREA[\"United States (USA) - onshore and offshore.\"], BBOX[15.56,167.65,74.71,-65.69]], ID[\"EPSG\",2163]]"
The default:
tf <- tempfile(fileext = '.tif')
(writeRaster(r, tf, overwrite=TRUE))
# coord. ref. : NAD27 / US National Atlas Equal Area (EPSG:9311)
Solution 1, usse the PROJ.4 notation
## get the proj.4 value
prj <- crs(r, proj=TRUE)
prj
#[1] "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +ellps=sphere +units=m +no_defs"
crs(r) <- prj
(writeRaster(r, tf, overwrite=TRUE))
#coord. ref. : +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +ellps=sphere +units=m +no_defs
Solution 2, use setGDALconfig. This requires terra version >= 1.5-27 (currently the development version, on windows you can install it with install.packages('terra', repos='https://rspatial.r-universe.dev'))
setGDALconfig("OSR_USE_NON_DEPRECATED", value="NO")
(writeRaster(r, tf, overwrite=TRUE))
# coord. ref. : US National Atlas Equal Area (EPSG:2163)