如何解决从PROJ4升级到PROJ6，并显示“废弃基准”警告 SF指向空间物体 SP指向空间物体光栅使用来自raster的自动重投影将点手动重新投影到栅格的crs 将栅格手动重新投影到点的crs

上下文

我的问题与从PROJ4升级到PROJ6引起的变化有关以及在各种R空间包（{{1}，sp，sf）中的后果。

我们现在收到许多有关“废弃原点”的警告，看起来有些令人担忧，我我对此应该怎么办有点疑惑。我看到这可能会带来可怕的后果在某些情况下，在其他情况下可以忽略它们。

似乎我不是唯一一个有点迷路的人（请参阅here）。我希望我的问题带有可复制的具体示例，有助于我们更好地理解该主题。

我了解我们可以remove the warnings，我已经阅读了上下文：r-spatial blog post， Migration to PROJ6/GDAL3 和these workshop notes（地图视图似乎在较新的版本中对此进行不同的处理）

问题

问题1：

可能是一个幼稚的问题：

我了解需要实施一种新的符号/格式（WKT）在PROJ6中（例如，由于需要更高的精度），但我不明白为什么会有需要从旧的proj4字符串符号中删除基准部分。为什么不只是保持原样（并以新的WKT格式/符号实现新功能）

问题2：

似乎有3种情况涉及旧proj4格式的数据丢失：

无警告：原点保留旧的proj4string表示法（rasterdéfault吗？）
我们有一个警告“已丢弃基准（...），但+ towgs84 =保留了值”
我们有一个警告“废弃的基准（...）”（在这种情况下，“ + towgs84 =” 字符串被丢弃/丢弃–> sf默认？）

下面的示例说明了我们收到这些警告的不同情况。

为什么在同一CRS（此处为epsg 31370）上有这3种不同的情况？
删除基准点和/或sp部分有什么后果？
我应该对第二个警告比对第三个警告更担心吗？

问题3：

在下面的可复制示例中，我尝试从对应的栅格中提取值到3点，栅格和点具有不同的CRS。但是，根据使用的方法，我会得到不同的结果。我的印象是，这与PROJ4 –> PROJ6的更改和基准有关下降，但我可能错了。我创建此示例仅是因为我想了解这些“基准下降”的后果在crs

我使用函数+towgs84和3种不同的通用方法（每次都使用 raster::extract和sf个对象），我希望从中得到相同的输出：

无需人工重新投影，我让sp进行工作以匹配点和栅格的crs
手动将点投影到栅格的crs
将栅格手动投影到点的crs

使用第三种方法，我为raster::extract或sp对象获得了2组不同的值，并且我用第二种方法（也是第一种方法）得到了第三组值（然后，如果我使用sf或sp对象）。

哪个值是“正确的”值？（可能是sf）
为什么某些方法会失败（这与警告消息/基准数据有关吗？）？
当我对空间物体进行投影并且收到这些警告时，如何检查我的空间物体已经“正确”投影了？

问题4：

是否可以就应该做什么提供一般性建议当我们收到这些警告消息时？

例如：

如果可能的话，不要使用旧的proj4字符串表示法，而应使用WKT表示法（但这并不总是可能的。例如，在这里-如果我没记错，我被迫使用旧的proj4表示法，因为这是351.7868 236.4216 309.0073的用法）
如果我知道这里使用的栅格的epsg 代码，我想最安全的方法是使用带有raster的代码重新投影点，例如：sf

可复制的示例

创建点和栅格空间数据+检查CRS的处理方式

SF指向空间物体

简而言之：CRS主要以WKT格式存储。旧的proj4string可应要求提供，并且不会删除基准/ st_transform(SF,crs = xxxx)部分

towgs84

SP指向空间物体

简而言之：CRS主要以proj4字符串格式存储，并删除基准点和library(sf) #> Linking to GEOS 3.8.0,GDAL 3.0.2,PROJ 6.2.1 library(sp) library(raster) # create 3 points coo <- data.frame(x = c(246000,247000,246500),y = c(184000,186000,185000)) # create an sf spatial object SF <- st_as_sf(coo,coords = c("x","y"),crs = 31370) # Check the CRS : # the proj4string includes the datum/+towgs84 information - no warning st_crs(SF)$proj4string #> [1] "+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-99.059,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs" # Same value with `raster::crs` but with a # Warning "discarded datum (...) but +towgs84= values preserved" raster::crs(SF) #> Warning in showSRID(uprojargs,format = "PROJ",multiline = "NO"): discarded datum UnkNown_based_on_International_1909_Hayford_ellipsoid in CRS deFinition,#> but +towgs84= values preserved #> CRS arguments: #> +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 #> +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl #> +towgs84=-99.059,-1 +units=m +no_defs # WKT st_crs(SF) #> Coordinate Reference System: #> User input: epsg:31370 #> wkt: #> BOUNDCRS[ #> SOURCECRS[ #> PROJCRS["Belge 1972 / Belgian LAmbert 72",#> BASEGEOGCRS["Belge 1972",#> DATUM["Reseau National Belge 1972",#> ELLIPSOID["International 1924",6378388,297,#> LENGTHUNIT["metre",1]]],#> PRIMEM["Greenwich",#> ANGLEUNIT["degree",0.0174532925199433]],#> ID["epsg",4313]],#> CONVERSION["Belgian LAmbert 72",#> METHOD["LAmbert Conic Conformal (2SP)",#> #> (...) #> #> ID["epsg",1609],#> REMARK["Scale difference is given by information source as 0.999999. Given in this record in ppm to assist application usage. Very similar parameter values (to slightly less precision) used for BD72 to ETRS89: see code 1652."]]] cat(raster::wkt(SF)) # does not work with sf #> Error in raster::wkt(SF): tentative d'obtenir le slot "crs" d'un objet (classe "sf") qui n'est pas un objet S4部分与towgs84不同）。新的WKT表示法将作为注释存储在CRS对象中，但与sf

不同

sf

光栅

在SP <- coo coordinates(SP) <- ~x+y proj4string(SP) <- CRS("+init=epsg:31370") #> Warning in showSRID(uprojargs,multiline = "NO"): discarded #> datum Reseau_National_Belge_1972 in CRS deFinition # the proj4 string do not contain the `towgs84` part # Warning "discarded datum (...)" CRS("+init=epsg:31370") #> Warning in showSRID(uprojargs,multiline = "NO"): discarded #> datum Reseau_National_Belge_1972 in CRS deFinition #> CRS arguments: #> +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 #> +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m #> +no_defs # With `raster::crs` same proj4string but no warning raster::crs(SP) #> CRS arguments: #> +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 #> +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m #> +no_defs # WKT notation + a warning: this WKT is indeed different from the SF one (no datum here ?) cat(comment(CRS("+init=epsg:31370"))) #> Warning in showSRID(uprojargs,multiline = "NO"): discarded #> datum Reseau_National_Belge_1972 in CRS deFinition #> PROJCRS["Belge 1972 / Belgian LAmbert 72",#> BASEGEOGCRS["Belge 1972",#> DATUM["Reseau National Belge 1972",#> ELLIPSOID["International 1924",#> LENGTHUNIT["metre",#> #> (...) #> #> USAGE[ #> ScopE["unkNown"],#> AREA["Belgium - onshore"],#> BBox[49.5,2.5,51.51,6.4]]] # Same output without warning cat(raster::wkt(SP)) #> PROJCRS["Belge 1972 / Belgian LAmbert 72",6.4]]]中，似乎同时存储了旧的proj4表示法和新的WKT表示法。

raster

此数据集的proj4string中不包含r <- raster::raster(system.file("external/test.Grd",package="raster")) raster::crs(r) #> CRS arguments: #> +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 #> +k=0.9999079 +x_0=155000 +y_0=463000 +datum=wgs84 +units=m +no_defs cat(raster::wkt(r)) #> PROJCRS["unkNown",#> BASEGEOGCRS["unkNown",#> DATUM["World Geodetic System 1984",#> ELLIPSOID["WGS 84",6378137,298.257223563,1]],#> ID["epsg",6326]],#> #> (...) #> #> AXIS["(N)",north,#> ORDER[2],#> LENGTHUNIT["metre",1,9001]]]]部分。但是当你读在proj4string中带有+towgs84部分的栅格似乎已被删除。
不可复制的示例：

+towgs84

我可能还应该探讨当我们使用GISfolder <- "/my/path" tmp <- raster(paste0(GISfolder,'my_file.tif')) #> Warning in showSRID(uprojargs,multiline = "NO"): discarded #> datum UnkNown_based_on_International_1909_Hayford_ellipsoid in CRS deFinition raster::crs(tmp) #> CRS arguments: #> +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=49.8333339 #> +lat_2=51.1666672333333 +x_0=150000.01256 +y_0=5400088.4378 +ellps=intl #> +units=m +no_defs cat(raster::wkt(tmp)) #> PROJCRS["unkNown",#> DATUM["UnkNown based on International 1909 (Hayford) ellipsoid",#> ELLIPSOID["International 1909 (Hayford)",#> ID["epsg",9001]]]],#> PRIMEM["Greenwich",#> ANGLEUNIT["degree",0.0174532925199433],8901]]],9001]]]]软件包而不是 stars，但这个问题已经很长了（我在光栅包上有很多代码）

从栅格中提取值

使用来自`raster`的自动重投影

raster::extract

将点手动重新投影到栅格的crs

光栅中带有proj4string的SF

# extract the values from the raster,# the function extract reprojects the points
# in the same crs as the raster layer
extract(r,SF)
#> Warning in showSRID(uprojargs,#>  but +towgs84= values preserved
#> Warning in .local(x,y,...): Transforming SpatialPoints to the CRS of the
#> Raster
#> [1] 351.7868 236.4216 309.0073
extract(r,SP)
#> Warning in .local(x,...): Transforming SpatialPoints to the CRS of the
#> Raster
#> [1] 351.7868 236.4216 309.0073

栅格中的WKT进行SF

SF_proj <- st_transform(SF,crs = raster::crs(r))
extract(r,SF_proj)
#> [1] 351.7868 236.4216 309.0073

光栅中带有proj4string的SP

SF_proj <- st_transform(SF,crs = raster::wkt(r))
extract(r,SF_proj)
#> [1] 351.7868 236.4216 309.0073

栅格中带有WKT的SP

SP_proj <- spTransform(SP,raster::crs(r)) extract(r,SP_proj) #> [1] 351.7868 236.4216 309.0073不接受wat fromat –>不起作用

sp::spTransform

将栅格手动重新投影到点的crs

使用SF的proj4string（带基准）

–>结果与以前的尝试不同

# error
SP_proj <- spTransform(SP,raster::wkt(r))
#> Error in CRS(CRSobj): PROJ4 argument-value pairs must begin with +: PROJCRS["unkNown",9001]]]]
# extract(r,SP_proj)

使用SF的WKT

–>不起作用，因为# epsg 31370 proj4 string with the datum: lAmbert72 <- sf::st_crs(31370)$proj4string lAmbert72 #> [1] "+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-99.059,-1 +units=m +no_defs" # there is a warning when we project the raster but the full string seems to be used r2 <- raster::projectRaster(r,crs = lAmbert72) #> Warning in showSRID(uprojargs,#> but +towgs84= values preserved raster::crs(r2) #> CRS arguments: #> +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 #> +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl #> +towgs84=-99.059,-1 +units=m +no_defs extract(r2,SP) #> [1] 341.6399 222.1028 301.2286不接受WKT格式其raster::projectRaster参数

crs

使用SP的proj4string（无基准）

–>结果与之前的尝试不同

lAmbert72 <- sf::st_crs(31370)
lAmbert72
#> Coordinate Reference System:
#>   User input: epsg:31370 
#>   wkt:
#> BOUNDCRS[
#>     SOURCECRS[
#>         PROJCRS["Belge 1972 / Belgian LAmbert 72",#>
#> (...)
#> 
#>             AREA["Belgium - onshore"],#>             BBox[49.5,6.4]],#>         ID["epsg",#>         REMARK["Scale difference is given by information source as 0.999999. Given in this record in ppm to assist application usage. Very similar parameter values (to slightly less precision) used for BD72 to ETRS89: see code 1652."]]]
r2 <- raster::projectRaster(r,crs = lAmbert72)
#> Error in wkt(projto): tentative d'obtenir le slot "crs" d'un objet (classe "crs") qui n'est pas un objet S4

sessionInfo

# epsg 31370 proj4 string without the datum:
lAmbert72 <- sp::CRS("+init=epsg:31370")@projargs
#> Warning in showSRID(uprojargs,multiline = "NO"): discarded
#> datum Reseau_National_Belge_1972 in CRS deFinition
lAmbert72
#> [1] "+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m +no_defs"
# warning
r3 <- raster::projectRaster(r,multiline = "NO"): discarded
#> datum UnkNown_based_on_International_1909_Hayford_ellipsoid in CRS deFinition
raster::crs(r3)
#> CRS arguments:
#>  +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
#> +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m
#> +no_defs

extract(r3,coo)
#> [1] 348.5775 329.1199 277.2260

用于：

sessionInfo()
#> R version 3.6.2 (2019-12-12)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.4 LTS
#>
#> (...)
#> 
#> other attached packages:
#> [1] raster_3.3-13 sp_1.4-2      sf_0.9-5      knitr_1.29   
#> 
#> (...)
#>

^{由reprex package（v0.3.0）于2020-09-03创建}

解决方法

https://gis.stackexchange.com/questions/372692中给出的一些注释。请首先去那里。

我了解需要新的符号/格式（WKT）在PROJ6中实现（例如，由于需要更高的精度），但是我不明白为什么需要从旧的proj4字符串表示法。为什么不保留它一如既往的（并以新的WKT格式/符号实现新功能）

GDAL中的+datum=部分已从GDAL> = 3中弃用。由于 sf ， rgdal 和 raster 如果使用GDAL来读取文件，则除了WGS84，NAD83和NAD27之外，Proj4字符串表示可能没有全部exportToProj4()。警告来自检查运行+datum=之前内部存在哪些节点以及之后存在哪些节点。当我们使用PROJ> = 6 / GDAL> = 3时，我们不能依靠exportToProj4()和+datum=。

有关示例的其他评论

+towgs84=

我正在使用开发版本以及PROJ和GDAL的最新版本。

> library(sf)
Linking to GEOS 3.8.1,GDAL 3.1.3,PROJ 7.1.1
> #> Linking to GEOS 3.8.0,GDAL 3.0.2,PROJ 6.2.1
> library(sp)
> library(raster)
> packageVersion("sf")
[1] ‘0.9.6’
> packageVersion("sp")
[1] ‘1.4.4’
> packageVersion("raster")
[1] ‘3.3.13’
> library(rgdal)
rgdal: version: 1.5-17,(SVN revision 1060)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.1.3,released 2020/09/01
Path to GDAL shared files: /usr/local/share/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 7.1.1,September 1st,2020,[PJ_VERSION: 711]
Path to PROJ shared files: /home/rsb/.local/share/proj:/usr/local/share/proj:/usr/local/share/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.4-4
To mute warnings of possible GDAL/OSR exportToProj4() degradation,use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

现在Proj4字符串中没有> coo <- data.frame(x = c(246000,247000,246500),y = c(184000,186000,185000)) > SF <- st_as_sf(coo,coords = c("x","y"),crs = 31370) > st_crs(SF)$proj4string [1] "+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-99.059,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs" > st_crs(SF) Coordinate Reference System: User input: EPSG:31370 wkt: PROJCRS["Belge 1972 / Belgian Lambert 72",BASEGEOGCRS["Belge 1972",DATUM["Reseau National Belge 1972",ELLIPSOID["International 1924",6378388,297,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4313]],CONVERSION["Belgian Lambert 72",METHOD["Lambert Conic Conformal (2SP)",9802]],PARAMETER["Latitude of false origin",90,0.0174532925199433],8821]],PARAMETER["Longitude of false origin",4.36748666666667,8822]],PARAMETER["Latitude of 1st standard parallel",51.1666672333333,8823]],PARAMETER["Latitude of 2nd standard parallel",49.8333339,8824]],PARAMETER["Easting at false origin",150000.013,1],8826]],PARAMETER["Northing at false origin",5400088.438,8827]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],1]],AXIS["northing (Y)",north,ORDER[2],USAGE[ SCOPE["unknown"],AREA["Belgium - onshore"],BBOX[49.5,2.5,51.51,6.4]],31370]]，但是WKT2_2019字符串中存在所有CRS规范。 +datum=对象中没有$proj4string，如果需要，它是即时生成的。

我们仍在进行强制，但已经有了：

"crs"

下一步：

> cat(raster::wkt(as(SF,"Spatial")),"\n")
PROJCRS["Belge 1972 / Belgian Lambert 72",31370]]

您注意到> SP <- coo > coordinates(SP) <- ~x+y > proj4string(SP) <- CRS("+init=epsg:31370") Warning message: In showSRID(uprojargs,format = "PROJ",multiline = "NO",prefer_proj = prefer_proj) : Discarded datum Reseau_National_Belge_1972 in CRS definition > cat(wkt(SP),8827]],19961]],AXIS["(E)",1,9001]]],AXIS["(N)",6.4]]]已经消失了，这是因为WKT2_2019中的DATUM绝对足以在需要时生成坐标操作。 PROJ> = 6 / GDAL> = 3无需转换为WGS84 GEOGCRS集线器，也无需转换为目标CRS。发出警告是因为+towgs84=既生成WKT2_2019字符串（已完全指定），又生成旧版Proj4字符串-现代PROJ / GDAL缺少位，我们希望没人再依赖它-如果这样做，您已经被警告。

我现在将其留在这里，指的是SE线程上的回复。如果 raster 开发人员可以发表评论，这将是有帮助的，但是据我们从反向依赖性检查中可以看到，光栅似乎已经过渡为在使用Proj4时优先使用WKT2_2019（与其他软件包一样）。 PROJ> = 6 / GDAL> =3。由于某些平台仍然是PROJ

根据我现在的理解得出部分答案。
NB：我完全不完全确定这一点。因此，如果我错了，请提供反馈…

总体思路是，sf，sp倾向于默认使用新的WKT表示法（可正确处理基准），即使它们可以显示（或根据要求检索）旧的且已弃用的proj4字符串表示法（带或不带基准）。

到目前为止（至少对我来说）情况尚不清楚关于能够提供WKT标记（raster）作为{ 字符串，但似乎仍然严重依赖proj4字符串。

因此，在大多数情况下，除非您强制使用 proj4表示法。但是对于raster::wkt，我仍然感到困惑，可能想念一些东西…… 我暂时不愿意使用raster。

我们现在可以尝试了解哪些答案是正确的，以及为什么：

raster::projectRaster

以下方法似乎是安全的，因为我们从栅格（作为library(sf) #> Linking to GEOS 3.8.0,PROJ 6.2.1 library(sp) library(raster) # create a raster r <- raster::raster(system.file("external/test.grd",package="raster")) # create an sf spatial object coo <- data.frame(x = c(246000,185000)) SF <- st_as_sf(coo,crs = 31370) # create an equivalent sp object SP <- coo coordinates(SP) <- ~x+y proj4string(SP) <- CRS(SRS_string = "EPSG:31370") # better than CRS("+init=epsg:31370") ??提供的字符串），我们将 raster::wkt和sf指向这个新的坐标参考系统。

sp

当我们使用SF_to_r <- st_transform(SF,crs = raster::wkt(r)) raster::extract(r,SF_to_r) # result (correct) : 351.7868 236.4216 309.0073 # note the use of `SRS_string` argument. `CRS(raster::wkt(r))` won't work SP_to_r <- spTransform(SP,CRS(SRS_string = raster::wkt(r))) raster::extract(r,SF_to_r) # result (correct) : 351.7868 236.4216 309.0073 class(raster::wkt(r)) # character时，以下方法似乎也可以工作（相同的结果）从raster::crs包中返回CRS类对象。我想这是因为 sp和sf均安全使用此对象提供的新WKT表示法（即使如果对象显然只包含proj4字符串，则WKT有点“隐藏” 在附加到对象的注释中）

sp

当我们投影栅格时，我们当然可以期望2 以下方法（一种使用SF_to_r <- st_transform(SF,crs = raster::crs(r)) raster::extract(r,SF_to_r) # result : 351.7868 236.4216 309.0073 SP_to_r <- spTransform(SP,raster::crs(r)) raster::extract(r,SF_to_r) # result : 351.7868 236.4216 309.0073 class(raster::crs(r)) # CRS class from `sp` str(raster::crs(r)) # 1 slot with the proj4 string cat(comment(raster::crs(r))) # this is where the WKT notation is "hidden"，一种使用sp），因为我们强制使用 proj4字符串（带有提供简单字符向量的sf和$proj4string）。应该始终避免这种情况……

我不确定这两个选项为何提供不同的结果但是我现在非常有信心这两个结果都是错误的。也许他们有所不同，因为数据在管道中的不同时刻被删除（提供的初始字符串为 sp和sf之间的字符向量是否不同？

@projargs

我们可以期望提供完整的CRS对象（而不是强制使用字符proj4字符串）即可解决此问题。但事实并非如此。也许是因为r_to_sf <- raster::projectRaster(r,crs = sf::st_crs(31370)$proj4string) raster::extract(r_to_sf,SF) # result (wrong) : 341.6399 222.1028 301.2286 r_to_sp <- raster::projectRaster(r,crs = sp::CRS("+init=epsg:31370")@projargs) raster::extract(r_to_sp,SF) # result (wrong) : 348.5775 329.1199 277.2260 class(sf::st_crs(31370)$proj4string) # character class(sp::CRS("+init=epsg:31370")@projargs) # character内部依赖于旧的proj4字符串??

不过，据罗杰·比万德（Roger Bivand）说：

从反向依赖检查中可以看到，栅格似乎已经过渡 PROJ> = 6 / GDAL> = 3
时，优先使用WKT2_2019（与其他软件包一样）优先于Proj4

所以我可能在某个地方错了，我仍然不知道如何安全地重新投影栅格对象...

raster

使用r_to_sp <- raster::projectRaster(r,crs = sp::CRS("+init=epsg:31370")) raster::extract(r_to_sp,SP) # result (wrong) : 341.6399 222.1028 301.2286 # same result with a slightly different syntax for CRS r_to_sp <- raster::projectRaster(r,crs = sp::CRS(SRS_string = "EPSG:31370")) raster::extract(r_to_sp,SP) # result (wrong) : 341.6399 222.1028 301.2286包，我们可以通过重新投影点或将点重新投影来获得“正确”结果。栅格。但是，似乎stars函数具有一些不可用的功能直接raster::extract（例如，使用多边形时为每个单元格计算权重）

Usefull raster vs stars function comparison

stars

这可能对反射有用吗？

Recommenadations from Roger Bivand：

如果可能，请避免使用Proj4字符串实例化“ CRS”对象，而应使用library(stars) STARS <- stars::read_stars(system.file("external/test.grd",package="raster")) # reproject the points into the same crs as the stars raster SF_to_STARS <- st_transform(SF,crs = st_crs(STARS)) aggregate(STARS,SF_to_STARS,function(x) x[1],as_points = FALSE)$test.grd # result (correct) = 351.7868 236.4216 309.0073 # reproject the stars raster into the same crs as the points STARS_to_SF <- st_transform(STARS,crs = st_crs(SF)) aggregate(STARS_to_SF,SF,as_points = FALSE)$test.grd # result (correct) = 351.7868 236.4216 309.0073

例如：

CRS(SRS_string=

所以（??）也可能是

# preferd syntax :
CRS(SRS_string = "OGC:CRS84")
#> Error in if (in_format == 4L) {: valeur manquante là où TRUE / FALSE est requis
# instead of :
CRS("+proj=longlat +datum=WGS84")

而不是：

CRS(SRS_string = "EPSG:31370")

避免使用CRS("+init=epsg:31370")，而是选择proj4string(x) <- proj4string(y)

^{由reprex package（v0.3.0）于2020-09-09创建}

这是我对您问题1的答复：

我了解需要新的符号/格式（WKT）在PROJ6中实现（例如，由于需要更高的精度），但是我不明白为什么需要从旧的proj4字符串表示法。为什么不保留它一如既往的（并以新的WKT格式/符号实现新功能）

我不知道为什么PROJ开发人员决定取消向后兼容性，但我认为这样做有充分的理由。而且没有人自愿参加这项工作。

作为R /空间开发人员（以及其他使用PROJ构建软件的人员），我们必须忍受这一点。问题是我们需要适应不同的PROJ版本（尤其是在Linux系统上）。试图向后兼容的同时向前移动会造成严重的混乱。

在R之类的脚本环境中，无法使用proj4表示法是一个真正的损失。proj4表示法可以直接理解； EPSG代码不透明，使用它们很容易导致错误。另外，如果没有可用的EPSG代码，则需要弄清楚如何编写自己的WKT。

栅格中的CRS

raster中的CRS对象与sp和rgdal中的CRS对象相同。它存储proj4和wkt表示法。罗杰·比万德（Roger Bivand）解释了为什么给出警告。

提取

要从栅格中提取值，请始终变换点（线，多边形），而不是栅格。转换栅格将导致新估计的值与原始值不同。请参阅讨论here