raster
raster copied to clipboard
rspatial
projectExtent error after moving script from Win Desktop (terra 1.6-17) to Centos Cluster (terra 1.6-21) | input crs is not valid
My code runs on my windows machine but not the Centos cluster. Any ideas why? It seems that something changed from terra 1.6-17 to terra 1.6-21.
library(raster)
# Define and project extent coordinates for DEM and NWM
DEM <- raster(paste0(wd2, "DEM_processed_OnR_1000m_NAD83_UTM11N.tif")) # Adjusted for PF overland
DEMext_org <- extent(DEM)
NWM_prj <- "+proj=lcc +units=m +a=6370000.0 +b=6370000.0 +lat_1=30.0 +lat_2=60.0 +lat_0=40.0 +lon_0=-97.0 +x_0=0 +y_0=0 +k_0=1.0 +nadgrids=@null +wktext +no_defs"
NWMext_org <- extent(matrix(data = c(-2303499.250, 2303500.750, -1919500.375, 1919499.625), byrow = T, ncol = 2, nrow = 2))
DEMext_prj <- extent(projectExtent(object = DEM, crs = crs(NWM_prj)))
#Error in h(simpleError(msg, call)) :
# error in evaluating the argument 'x' in selecting a method for function 'extent': [project] input crs is not valid
Session on Centos:
> sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.3.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
[4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] retry_0.1.0 doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2 raster_3.6-3
[6] rgdal_1.5-32 sp_1.5-0 ncdf4_1.19 aws.s3_0.3.21
loaded via a namespace (and not attached):
[1] Rcpp_1.0.9 lattice_0.20-45 codetools_0.2-18 digest_0.6.29
[5] aws.signature_0.6.0 grid_4.2.1 R6_2.5.1 httr_1.4.4
[9] rlang_1.0.5 curl_4.3.2 xml2_1.3.3 tools_4.2.1
[13] compiler_4.2.1 terra_1.6-21 base64enc_0.1-3
Session on Windows:
> sessionInfo()
R version 4.2.1 (2022-06-23 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C LC_TIME=English_United States.utf8
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2 raster_3.6-3 rgdal_1.5-32 sp_1.5-0 tictoc_1.1
[8] retry_0.1.0 ncdf4_1.19 aws.s3_0.3.21
loaded via a namespace (and not attached):
[1] Rcpp_1.0.9 lattice_0.20-45 codetools_0.2-18 digest_0.6.29 aws.signature_0.6.0 grid_4.2.1
[7] R6_2.5.1 httr_1.4.4 cli_3.4.0 rlang_1.0.5 curl_4.3.2 xml2_1.3.3
[13] tools_4.2.1 compiler_4.2.1 terra_1.6-17 base64enc_0.1-3
It is probably related to the version of GDAL/PROJ . What does
terra::gdal(lib="")
return on CENTOS?
And this might work:
DEMext_prj <- extent(projectExtent(object = DEM, crs = NWM_prj)
(or using terra instead of raster)
@rhijmans, these are the GDAL/PROJ versions: on Centos
> rgdal::getGDALVersionInfo()
[1] "GDAL 2.3.2, released 2018/09/21"
> rgdal::getPROJ4VersionInfo()
[1] "Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]"
attr(,"short")
[1] 480
extent(projectExtent(object = DEM, crs = NWM_prj)) did not work:
Error in h(simpleError(msg, call)) :
error in evaluating the argument 'x' in selecting a method for function 'extent': [project] input crs is not valid```
These are very old versions. But I can run this with nearly as old versions of GDAL/PROJ on Debian.
library(raster)
DEM <- raster()
NWM_prj <- "+proj=lcc +units=m +a=6370000.0 +b=6370000.0 +lat_1=30.0 +lat_2=60.0 +lat_0=40.0 +lon_0=-97.0 +x_0=0 +y_0=0 +k_0=1.0 +nadgrids=@null +wktext +no_defs"
extent(projectExtent(object = DEM, crs = crs(NWM_prj)))
#class : Extent
#xmin : -292426517
#xmax : 88725532
#ymin : -166378680
#ymax : 60166087
terra::gdal(lib="")
# gdal proj geos
#"2.4.0" "5.2.0" "3.7.1"
packageVersion("raster")
#[1] ‘3.6.3’
packageVersion("terra")
#[1] ‘1.6.17’
I assume this also fails for you?
It failed too. Here are the versions:
> terra::gdal(lib="")
gdal proj geos
"2.3.2" "4.8.0" "3.5.0"
> packageVersion("raster")
[1] ‘3.6.3’
> packageVersion("terra")
[1] ‘1.6.21’
Here is the traceback:
> DEMext_prj <- extent(projectExtent(object = DEM, crs = crs(NWM_prj)))
Error in h(simpleError(msg, call)) :
error in evaluating the argument 'x' in selecting a method for function 'extent': [project] input crs is not valid
> traceback()
11: h(simpleError(msg, call))
10: .handleSimpleError(function (cond)
.Internal(C_tryCatchHelper(addr, 1L, cond)), "[project] input crs is not valid",
base::quote(NULL))
9: stop("[", f, "] ", emsg, ..., call. = FALSE)
8: error(f, x@ptr$getError())
7: messages(x, "project")
6: .local(x, ...)
5: terra::project(xy, projto)
4: terra::project(xy, projto)
3: .rawTransform(projfrom, projto, xy)
2: projectExtent(object = DEM, crs = crs(NWM_prj))
1: extent(projectExtent(object = DEM, crs = crs(NWM_prj)))
Thanks!
Do the examples in ?projectRaster work? To see if this is a general failure (I assume it is) or is it associated with your specific crs).
Can you get more up to date GDAL/PROJ libraries installed ?
@rhijmans, I already emailed the administrator to see if GDAL/PROJ can be updated. If not, do you know if it is possible to install GDAL/PROJ in my home directory and make R/RStudio identify them?
The examples in ?projectRaster work perfectly. I even ran the projectRaster example using my crs and it worked. Do you still think the problem is the crs due to the versions of GDAL/PROJ? I think there is something odd with projectExtent.
I just realized that the problem is that projectExtent cannot extract the crs of my DEM:
> crs(DEM)
Coordinate Reference System:
Deprecated Proj.4 representation: NA
Warning message:
In wkt(x) : CRS object has no comment
While on Windows:
> crs(DEM)
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
WKT2 2019 representation:
PROJCRS["NAD83 / UTM zone 11N",
BASEGEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
CONVERSION["UTM zone 11N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-117,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["North America - between 120°W and 114°W - onshore and offshore. Canada - Alberta; British Columbia; Northwest Territories; Nunavut. United States (USA) - California; Idaho; Nevada, Oregon; Washington."],
BBOX[30.88,-120,83.5,-114]],
ID["EPSG",26911]]
Thanks, that is helpful, that explains it. But why? What does this return?
proj4string(DEM)
And this:
proj4string(DEM) <- "+proj=utm +zone=11 +datum=NAD83"
proj4string(DEM)
You cannot install your own GDAL/PROJ libs, but you could install an older version (the previous release) of "raster" that uses rgdal instead of terra for raster handling.
Here the results:
> proj4string(DEM)
[1] NA
> proj4string(DEM) <- "+proj=utm +zone=11 +datum=NAD83"
> proj4string(DEM)
[1] "+proj=utm +zone=11 +datum=NAD83 +ellps=GRS80 +towgs84=0,0,0"
Which raster version should I install and how? Thanks!
See here for installing older versions of R packages. The versions of "raster" are listed here. Version 3.5-21 should be safe (not affected by the changes from moving away from rgdal)
Your results suggest that if you do
proj4string(DEM) <- "+proj=utm +zone=11 +datum=NAD83"
your script may work again with the versions you have now.
proj4string(DEM) <- "+proj=utm +zone=11 +datum=NAD83" works for now. Thanks!