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

编程问答 2022-08-12

问题描述

上下文

我的问题与从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格式的数据丢失：

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

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

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

问题3：

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

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

  1. 无需人工重新投影，我让sp进行工作以匹配点和栅格的crs
  2. 手动将点投影到栅格的crs
  3. 将栅格手动投影到点的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 @R_176_4045@ion - 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 @R_176_4045@ion 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 @R_176_4045@ion 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中给出的一些注释。请首先去那里。

  1. 我了解需要新的符号/格式（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

rrrasterrgdalsfsp