我正在尝试使用terra包读取几个光栅;首先,我得到一个带有文件名的向量,然后我使用rast()
函数读取文件。
Land_variables_files <- list.files("Land_Variables1", pattern = '.tif$', full.names = T)
Land_variables_temp <- rast(Land_variables_files)
它们应该已经在EPSG 3035中了事实上就是这样,但是当我写它们的时候我在屏幕上看到了这个信息
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
我也试着在EPSG 4326中阅读它们,但没有任何变化。
可能是什么问题?我的参考系统或我设置它的方式有错误?我在QGIS中进行初步分析,并将光栅保存在EPSG 3035中。
使用show(Land_Variables_temp)
,这是我得到的:
class : SpatRaster
dimensions : 1307, 1307, 7 (nrow, ncol, nlyr)
resolution : 100, 100 (x, y)
extent : 4496516, 4627216, 2494261, 2624961 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
sources : Agricultural_Fields.tif
Coniferous_Forests.tif
Deciduous_Forests.tif
... and 4 more source(s)
names : Agric~ields, Conif~rests, Decid~rests, Global_Model, Pastures, Settlements, ...
min values : ? , ? , ? , 64, ? , ? , ...
max values : ? , ? , ? , 100, ? , ? , ... '
我不知道为什么有些变量有?作为最小/最大值,但绘制它们或单独取它们的值是可以的。
另一方面,使用crs()
,这是结果:
[1] "PROJCRS["ETRS89-extended / LAEA Europe",n BASEGEOGCRS["ETRS89",n ENSEMBLE["European Terrestrial Reference System 1989 ensemble",n MEMBER["European Terrestrial Reference Frame 1989"],n MEMBER["European Terrestrial Reference Frame 1990"],n MEMBER["European Terrestrial Reference Frame 1991"],n MEMBER["European Terrestrial Reference Frame 1992"],n MEMBER["European Terrestrial Reference Frame 1993"],n MEMBER["European Terrestrial Reference Frame 1994"],n MEMBER["European Terrestrial Reference Frame 1996"],n MEMBER["European Terrestrial Reference Frame 1997"],n MEMBER["European Terrestrial Reference Frame 2000"],n MEMBER["European Terrestrial Reference Frame 2005"],n MEMBER["European Terrestrial Reference Frame 2014"],n ELLIPSOID["GRS 1980",6378137,298.257222101,n LENGTHUNIT["metre",1]],n ENSEMBLEACCURACY[0.1]],n PRIMEM["Greenwich",0,n ANGLEUNIT["degree",0.0174532925199433]],n ID["EPSG",4258]],n CONVERSION["Europe Equal Area 2001",n METHOD["Lambert Azimuthal Equal Area",n ID["EPSG",9820]],n PARAMETER["Latitude of natural origin",52,n ANGLEUNIT["degree",0.0174532925199433],n ID["EPSG",8801]],n PARAMETER["Longitude of natural origin",10,n ANGLEUNIT["degree",0.0174532925199433],n ID["EPSG",8802]],n PARAMETER["False easting",4321000,n LENGTHUNIT["metre",1],n ID["EPSG",8806]],n PARAMETER["False northing",3210000,n LENGTHUNIT["metre",1],n ID["EPSG",8807]]],n CS[Cartesian,2],n AXIS["northing (Y)",north,n ORDER[1],n LENGTHUNIT["metre",1]],n AXIS["easting (X)",east,n ORDER[2],n LENGTHUNIT["metre",1]],n USAGE[n SCOPE["Statistical analysis."],n AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Turkey; United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],n BBOX[24.6,-35.58,84.73,44.83]],n ID["EPSG",3035]]"
使用gdal(lib = "all")
作为回报,我得到这个:
gdal proj geos
"3.5.2" "8.2.1" "3.9.3"
提前感谢大家
它适合我
library(terra)
#terra 1.6.53
r <- rast(nrow=2, ncol=2, crs="EPSG:3035", vals=1:4)
x <- writeRaster(r, "test.tif", overwrite=T)
x
#class : SpatRaster
#dimensions : 2, 2, 1 (nrow, ncol, nlyr)
#resolution : 180, 90 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
#source : test.tif
#name : lyr.1
#min value : 1
#max value : 4
也许是因为我有不同版本的GDAL/PROJ库?我在这个例子中使用的是windows。
gdal(lib="all")
# gdal proj geos
# "3.6.2" "9.1.1" "3.11.1"
这个问题以前已经报告过了。到目前为止,似乎只有使用ETRS数据的crs才会给出此消息。这个问题没有解决,但似乎并没有产生什么问题。