logo标签列表

PROJ4升级到PROJ6和“废弃数据”警告

rrasterr-sfr-sprgdal
20

Context

这些问题与从PROJ4升级到PROJ6所引起的更改及其对各种R空间包（sp，sf，raster）的影响有关。

我们现在收到了很多关于“已放弃基准”的警告，看起来有点令人担忧，我有点困惑我应该做什么。我可以看到这可能会在某些情况下产生可怕的后果，并且在其他情况下可以忽略它们。

看来我不是唯一一个有点迷失的人（请参见此处）。我希望我的问题和具有特定可再现示例的实例将帮助我们更好地理解这个主题。

我理解我们可以消除警告，并且我已经阅读了上下文: r-spatial博客文章、 升级到PROJ6/GDAL3 和这些研讨会笔记（mapview似乎以不同的方式处理此问题在最近的版本中）。

Questions

问题1:

可能是一个天真的问题：

我理解需要一种新的符号/格式（WKT）在PROJ6中实现（例如由于需要更高的精度），但我不明白为什么需要将旧proj4字符串符号中的基准部分删除。为什么不保持它一直如此（并在新的WKT格式/符号中实现新功能）？

问题2:

似乎我们有3种情况涉及旧proj4格式中的基准删除：

  1. 没有警告：基准与旧proj4string符号相同（sf默认值？）
  2. 我们收到了一个警告“已放弃基准（…），但保留了+ towgs84 =值”
  3. 我们收到了一个警告“已放弃基准（…）”（在这种情况下，“+towgs84 =”字符串的一部分被丢弃/删除->sp默认值？）

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

为什么在相同的CRS（这里是EPSG 31370）上有这3种不同的情况？
删除基准和/或+towgs84部分的后果是什么？
第二个警告比第三个警告更少使我担心吗？

问题3:

在下面的可再现示例中，我尝试从栅格提取与3个点对应的值，栅格和点具有不同的CRS。然而，根据所使用的方法不同，我会得到不同的结果。我有一种印象，这与这个PROJ4 -> PROJ6更改和基准删除有关，但我可能错了。我创建

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,53.322,-112.486,0.419,-0.83,1.885,-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",0,
#>                     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

SP点空间对象

简要说，CRS主要以proj4字符串格式存储，并且不像sf那样保留数据和towgs84部分。新的WKT符号被存储为CRS对象中的注释，但与sf不同。

SP <- coo
coordinates(SP) <- ~x+y
proj4string(SP) <- CRS("+init=epsg:31370")
#> Warning in showSRID(uprojargs, format = "PROJ", 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, format = "PROJ", 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, format = "PROJ", 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",6378388,297,
#>                 LENGTHUNIT["metre",1]]],
#>
#> (...)
#> 
#>     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",
#>     BASEGEOGCRS["Belge 1972",
#>         DATUM["Reseau National Belge 1972",
#>             ELLIPSOID["International 1924",6378388,297,
#>                 LENGTHUNIT["metre",1]]],
#>
#> (...)
#> 
#>     USAGE[
#>         SCOPE["unknown"],
#>         AREA["Belgium - onshore"],
#>         BBOX[49.5,2.5,51.51,6.4]]]

栅格

在栅格中，既存的proj4表示法和新的WKT表示法似乎都被存储了。

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,
#>                 LENGTHUNIT["metre",1]],
#>             ID["EPSG",6326]],
#>
#> (...)
#> 
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]]]

这个数据集的proj4string部分中不包含+towgs84部分。但是当您读取具有proj4string中+towgs84部分的栅格时，似乎会被删除。
无法重现的示例：

GISfolder <- "/my/path"
tmp <- raster(paste0(GISfolder, 'my_file.tif'))
#> Warning in showSRID(uprojargs, format = "PROJ", 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",
#>     BASEGEOGCRS["unknown",
#>         DATUM["Unknown based on International 1909 (Hayford) ellipsoid",
#>             ELLIPSOID["International 1909 (Hayford)",6378388,297,
#>                 LENGTHUNIT["metre",1,
#>                     ID["EPSG",9001]]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8901]]],
#>
#> (...)
#> 
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]]]

我应该还需要探索使用stars包而不是raster时会发生什么，但这个问题已经很长了（而且我已经有很多基于raster包的代码）

从栅格中提取值

使用raster::extract的自动投影

# 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, format = "PROJ", multiline = "NO"): Discarded datum Unknown_based_on_International_1909_Hayford_ellipsoid in CRS definition,
#>  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, y, ...): Transforming SpatialPoints to the CRS of the
#> Raster
#> [1] 351.7868 236.4216 309.0073

手动将点重新投影到栅格数据的坐标参考系统

使用栅格数据的proj4string进行投影。
SF_proj <- st_transform(SF, crs = raster::crs(r))
extract(r, SF_proj)
#> [1] 351.7868 236.4216 309.0073
使用WKT从栅格中提取SF。
SF_proj <- st_transform(SF, crs = raster::wkt(r))
extract(r, SF_proj)
#> [1] 351.7868 236.4216 309.0073
从栅格中获取带有proj4string的SP对象。
SP_proj <- spTransform(SP, raster::crs(r))
extract(r, SP_proj)
#> [1] 351.7868 236.4216 309.0073

在栅格图像中使用WKT格式的空间参考信息，但是sp::spTransform不支持该格式，因此无法使用。

# error
SP_proj <- spTransform(SP, raster::wkt(r))
#> Error in CRS(CRSobj): PROJ4 argument-value pairs must begin with +: PROJCRS["unknown",
#>     BASEGEOGCRS["unknown",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]],
#>             ID["EPSG",6326]],
#>
#> (...)
#> 
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]]]
# extract(r, SP_proj)

手动将栅格转换到点的坐标参考系统

使用SF的proj4string（带有数据），

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

# 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,53.322,-112.486,0.419,-0.83,1.885,-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, format = "PROJ", multiline = "NO"): Discarded datum Unknown_based_on_International_1909_Hayford_ellipsoid in CRS definition,
#>  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,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs

extract(r2, SP)
#> [1] 341.6399 222.1028 301.2286
使用SF的WKT格式不起作用，因为raster::projectRaster不接受WKT格式作为其crs参数。
lambert72 <- sf::st_crs(31370)
lambert72
#> 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",0,
#>
#> (...)
#> 
#>             AREA["Belgium - onshore"],
#>             BBOX[49.5,2.5,51.51,6.4]],
#>         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."]]]
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
使用SP的proj4string（不带数据）会导致与之前尝试的结果不同。
# EPSG 31370 proj4 string without the datum:
lambert72 <- sp::CRS("+init=epsg:31370")@projargs
#> Warning in showSRID(uprojargs, format = "PROJ", 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, crs = lambert72)
#> Warning in showSRID(uprojargs, format = "PROJ", 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   
#> 
#> (...)
#>

用于：

GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1

由reprex包（v0.3.0）于2020年09月03日创建

3个回答
7

以下是https://gis.stackexchange.com/questions/372692中提到的一些注意事项，请先查看。

  1. 我知道需要新的符号/格式(WKT)在PROJ6中实现(例如，由于需要更高的精度)，但我不明白为什么要从旧proj4字符串符号中移除数据部分。为什么不像以前一样保留它(并在新的WKT格式/符号中实现新功能)?

+datum=部分在GDAL的exportToProj4()自GDAL >= 3中已经被弃用。由于sf、rgdal和raster使用GDAL读取文件，因此Proj4字符串表示法没有所有+datum=，也许除了WGS84、NAD83和NAD27。这些警告来自于在运行exportToProj4()之前和之后检查内部存在哪些节点。当我们使用PROJ>=6/GDAL>=3时，我们不能依赖于+ datum=和+towgs84=。

下面是进一步的评论:

> 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.

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

> 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",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4313]],
    CONVERSION["Belgian Lambert 72",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",90,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",4.36748666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",51.1666672333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",49.8333339,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",150000.013,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",5400088.438,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Belgium - onshore"],
        BBOX[49.5,2.5,51.51,6.4]],
    ID["EPSG",31370]]

现在Proj4字符串中不再存在+datum=，但所有CRS规范都包含在WKT2_2019字符串中。在"crs"对象中没有$proj4string，如果你需要，它会被即时生成。

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

> cat(raster::wkt(as(SF, "Spatial")), "\n")
PROJCRS["Belge 1972 / Belgian Lambert 72",
    BASEGEOGCRS["Belge 1972",
        DATUM["Reseau National Belge 1972",
            ELLIPSOID["International 1924",6378388,297,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4313]],
    CONVERSION["Belgian Lambert 72",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",90,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",4.36748666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",51.1666672333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",49.8333339,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",150000.013,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",5400088.438,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Belgium - onshore"],
        BBOX[49.5,2.5,51.51,6.4]],
    ID["EPSG",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), "\n")
PROJCRS["Belge 1972 / Belgian Lambert 72",
    BASEGEOGCRS["Belge 1972",
        DATUM["Reseau National Belge 1972",
            ELLIPSOID["International 1924",6378388,297,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4313]],
    CONVERSION["Belgian Lambert 72",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",90,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",4.36748666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",51.1666672333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",49.8333339,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",150000.013,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",5400088.438,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]],
        ID["EPSG",19961]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
    USAGE[
        SCOPE["unknown"],
        AREA["Belgium - onshore"],
        BBOX[49.5,2.5,51.51,6.4]]]

您注意到+towgs84 = 已经删除，这是因为在WKT2_2019中的DATUM足以生成需要的坐标操作。PROJ >= 6/GDAL >= 3无需转换为WGS84 GEOGCRS中心并继续转换到目标CRS。警告出现是因为 sp::CRS() 生成了完全指定的WKT2_2019字符串和传统的Proj4字符串——对于现代PROJ/GDAL缺少一些位——我们希望没有人再依赖它们——如果您这样做了，那么您已经被警告了。
我现在将这里留下来，参考SE线程上的回复。如果一个raster开发者能够评论，这将是有帮助的，但据我们所知，从反向依赖检查中可以看出，raster似乎已经过渡到在PROJ >= 6/GDAL >= 3情况下优先使用WKT2_2019（像其他软件包一样）而不是Proj4。由于某些平台仍然是PROJ < 6/GDAL < 3，因此我们必须尽可能提供两种设置。
谢谢您的回答。我现在对情况有了更好的理解。如果我没错的话，我的例子中的问题来自于raster::projectRaster函数，无论您提供的crs格式是什么（作为CRS sp对象或当然作为proj4字符串），它都无法正确投影栅格。我尝试在单独的答案中澄清每个选项的工作原理或失败原因（如果我错了，请纠正我）。 - Gilles San Martin
3
基于我现在理解的部分答案。
注意：我并不完全确定。如果我错了，请提供反馈...
总体想法是，sf、sp 默认使用新的 WKT 标记法（正确处理数据），即使可以显示（或根据请求检索）旧的和已弃用的 proj4 字符串标记法（带或不带数据）。
至少对我来说，关于raster的情况还不太清楚。它能够提供WKT标记法（raster::wkt）作为字符字符串，但似乎仍然严重依赖proj4字符串。
因此，大多数情况下，投影应该是正常的，除非您强制使用proj4标记法。但是对于raster，我仍然感到困惑，可能我还缺少一些东西...目前我会非常不舒服使用raster::projectRaster。
我们现在可以尝试理解哪些答案是正确的以及为什么：
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.2, 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, 247000, 246500), y = c(184000, 186000, 185000))
SF <- st_as_sf(coo, coords = c("x", "y"), 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") ??

以下方法似乎是安全的，因为我们从栅格中提取WKT符号（作为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时似乎也有效（结果相同），它返回一个sp包中的CRS类对象。我猜这是因为sf和sp都安全地使用了此对象可用的新WKT表示法（即使对象仅包含proj4字符串，WKT也被“隐藏”在附加到对象的注释中）。
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，另一种使用sf）会失败，因为我们强制使用proj4字符串（使用$proj4string和@projargs提供一个简单的字符向量）。这应该始终被避免...
我不确定为什么这两个选项会提供不同的结果，但我相当有信心现在两个结果都是错误的。也许它们之间的差异是因为在管道中的不同时刻丢弃了数据（最初提供的字符串在sp和sf之间作为字符向量是不同的）？
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

我们可以期望提供完整的CRS对象（而不是强制使用字符proj4字符串）来解决问题。但这似乎并不是情况。也许是因为raster在内部依赖于旧的proj4字符串？？
然而，根据Roger Bivand的说法：
据我们从反向依赖检查中看到的，raster似乎已经过渡到使用WKT2_2019（就像其他软件包一样），在PROJ>=6/GDAL>=3时优先于Proj4。
所以我可能在某个地方错了，我仍然不知道如何安全地重新投影光栅对象...
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 函数具有一些不容易直接使用 stars 的特性（例如，在使用多边形时计算每个单元格的权重）。

有用的 raster vs stars 函数比较 

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, function(x) x[1], as_points = FALSE)$test.grd # result (correct) = 351.7868 236.4216 309.0073
这可能对反射也有用吗？ Roger Bivand的建议：

如果可能的话，避免使用Proj4字符串实例化“CRS”对象，而是使用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)，而是优先使用： slot(x, "proj4string") <- slot(y, "proj4string")

此内容创建于2020年9月9日，由 reprex package (v0.3.0) 创建。

2
这是我对你的问题1的回答：
我知道在PROJ6中实现新的符号/格式（例如由于需要更高的精度），需要采用新的WKT符号/格式，但我不明白为什么需要从旧的proj4字符串符号中删除数据部分。为什么不像以前一样保留它（并在新的WKT格式/符号中实现新的功能）？
我不知道PROJ开发人员为什么决定打破向后兼容性，但我认为他们有非常充分的理由；而且没有人自愿去解决这个问题。
作为R / spatial开发人员（以及其他使用PROJ构建软件的人），我们必须应对这种情况。问题是我们需要适应不同的PROJ版本（特别是在Linux系统上）。在保持向后兼容的同时推进难以取得进展，造成了可怕的混乱。
在像R这样的脚本环境中无法使用proj4符号确实是一个真正的损失。proj4符号可以直接理解；EPSG代码是不透明的，容易导致错误。此外，如果没有EPSG代码可用，您需要找出如何编写自己的WKT。
栅格中的CRS

raster 中的 CRS 对象与 sp 和 rgdal 中的相同。它存储了 proj4 和 wkt 表示法。Roger Bivand 解释了为什么会出现警告。

提取

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

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接