那些遇到过的问题

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

Gil*_*les 13 r raster rgdal r-sp r-sf

语境

我的问题与从 PROJ4 升级到 PROJ6 引起的变化以及各种 R 空间包 ( sp, sf, raster) 中的后果有关。

我们现在收到很多关于“废弃数据”的警告，看起来有点令人担忧，我有点困惑我应该怎么做。我可以看到这在某些情况下会产生可怕的后果，而在其他情况下可以忽略它们。

似乎我不是唯一一个有点迷茫的人（见这里）。我希望我提出的带有特定可重现示例的问题将有助于我们更好地理解该主题。

我知道我们可以删除警告，并且我已经阅读了上下文：r-spatial 博客文章、 迁移到 PROJ6/GDAL3 和这些研讨会笔记（mapview 在更新的版本中似乎以不同的方式处理此问题）

问题

问题 1：

可能是一个幼稚的问题：

我知道需要在 PROJ6 中实现新的符号/格式（WKT）（例如，因为需要更高的精度），但我不明白为什么需要从旧的 proj4 字符串中删除数据部分符号。为什么不保持原样（并以新的 WKT 格式/符号实现新功能）

问题2 ：

似乎我们有 3 个关于旧 proj4 格式的数据丢失的案例：

  1. 没有警告：数据保持在旧的 proj4string 表示法中（sf默认？）
  2. 我们有一个警告“丢弃的数据 (...) 但 +towgs84= 保留了值”
  3. 我们有一个警告“Discarded datum (...)”（在这种情况下，字符串的“+towgs84=”部分被丢弃/丢弃 –>sp默认 ?? ）

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

为什么我们在同一个 CRS 上有这 3 个不同的案例（这里是 EPSG 31370）？
删除基准和/或+towgs84零件的后果是什么？
我应该对第二个警告比第三个更不担心吗？

问题 3：

在下面的可重现示例中，我尝试从对应于 3 个点的栅格中提取值，栅格和具有不同 CRS 的点。然而，根据使用的方法，我得到不同的结果。我的印象是这与 PROJ4 –> PROJ6 更改和数据丢失有关，但我可能是错的。我创建这个例子只是因为我想了解这些“数据丢失”在 crs 中的后果

我使用函数raster::extract和 3 种不同的一般方法（每次都 用于点sf和sp对象），我希望从中得到相同的输出：

  1. 没有手动重新投影，我让raster::extract工作来匹配点和栅格的 crs
  2. 手动将点重新投影到栅格的 crs
  3. 手动将栅格重新投影到点的 crs

使用第 3 种方法，我得到 2 组不同的sp或sf对象值，并且我得到了第二种（和第一种方法）的第三组值（如果我使用sporsf对象，则这些值是相同的）。

  • 哪些值是“正确的”值？（大概351.7868 236.4216 309.0073）
  • 为什么有些方法会失败（这与警告消息/数据内容有关吗？）？
  • 当我对空间对象进行投影并收到这些警告时，如何检查我的空间对象是否已“正确”投影？

问题 4：

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

例如 ：

  • 如果可能的话，不要使用旧的 proj4 字符串表示法，而是使用 WKT 表示法（但这并不总是可能的。例如 - 如果我没有错，我被迫使用旧的 proj4 表示法，因为这是raster使用的）
  • 如果我知道这里使用的栅格的 epsg 代码，我想最安全的方法是使用该代码重新投影点sf，例如：st_transform(SF, crs = xxxx)

可重现的例子

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

SF点空间对象

简而言之：CRS主要以WKT格式存储。旧的 proj4string 可应要求提供，不要丢弃数据/towgs84部分

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
Run Code Online (Sandbox Code Playgroud)

SP点空间对象

简而言之：CRS 主要以 proj4 字符串格式存储，并删除数据和towgs84部分，不像sf）。新的 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]]]
Run Code Online (Sandbox Code Playgroud)

光栅

在 中raster，似乎存储了旧的 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]]]]
Run Code Online (Sandbox Code Playgroud)

此数据集不包含+towgs84proj4string 中的部分。但是，当您读取+towgs84proj4string 中有一部分的光栅时，它似乎被删除了。
不可复制的例子：

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]]]]
Run Code Online (Sandbox Code Playgroud)

我可能还应该探索当我们使用stars包而不是时会发生什么， 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
Run Code Online (Sandbox Code Playgroud)

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

带有来自栅格的 proj4string 的 SF 
SF_proj <- st_transform(SF, crs = raster::crs(r))
extract(r, SF_proj)
#> [1] 351.7868 236.4216 309.0073
Run Code Online (Sandbox Code Playgroud) 来自栅格的带有 WKT 的 SF 
SF_proj <- st_transform(SF, crs = raster::wkt(r))
extract(r, SF_proj)
#> [1] 351.7868 236.4216 309.0073
Run Code Online (Sandbox Code Playgroud) 来自光栅的带有 proj4string 的 SP 
SP_proj <- spTransform(SP, raster::crs(r))
extract(r, SP_proj)
#> [1] 351.7868 236.4216 309.0073
Run Code Online (Sandbox Code Playgroud) 来自栅格的带有 WKT 的 SP

wkt fromat 不被接受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)
Run Code Online (Sandbox Code Playgroud)

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

使用 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
Run Code Online (Sandbox Code Playgroud) 使用SF的WKT

–> 不起作用，因为raster::projectRaster它的crs参数不接受 WKT 格式

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
Run Code Online (Sandbox Code Playgroud) 使用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
Run Code Online (Sandbox Code Playgroud)

会话信息

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   
#> 
#> (...)
#>
Run Code Online (Sandbox Code Playgroud)

用于：

GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
Run Code Online (Sandbox Code Playgroud)

由reprex 包(v0.3.0)于 2020 年 9 月 3 日创建

小智 6

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

  1. 我知道需要在 PROJ6 中实现新的符号/格式（WKT）（例如，因为需要更高的精度），但我不明白为什么需要从旧的 proj4 字符串中删除数据部分符号。为什么不保持原样（并以新的 WKT 格式/符号实现新功能）

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

进一步的评论与示例有关：

> 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.
Run Code Online (Sandbox Code Playgroud)

我正在使用开发版本，以及最新版本的 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]]
Run Code Online (Sandbox Code Playgroud)

现在+datum=Proj4 字符串中没有剩余，但所有 CRS 规范都存在于 WKT2_2019 字符串中。对象中没有$proj4string，"crs"如果您要求，它会即时生成。

我们仍在研究强制转换，但我们已经有了：

> 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]]
Run Code Online (Sandbox Code Playgroud)

下一个：


> 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]]]
Run Code Online (Sandbox Code Playgroud)

你注意到+towgs84=了，那是因为 WKT2_2019 中的 DATUM 绝对足以在需要时生成坐标操作。PROJ >= 6/GDAL >= 3 不需要转换到 WGS84 GEOGCRS hub 并向前转换到目标 CRS。发出警告是因为sp::CRS()生成了完全指定的 WKT2_2019 字符串和遗留的 Proj4 字符串——现代 PROJ/GDAL 缺少位，我们希望没有人再依赖它——如果你这样做了，你已经被警告了。

我现在把这个留在这里，参考 SE 线程上的回复。如果光栅开发人员可以发表评论，这会很有帮助，但就我们从反向依赖性检查中看到的而言，当 PROJ >= 6/GDAL > = 3.由于有些平台还是PROJ < 6/GDAL < 3，所以我们要尽量提供这两种设置。

归档时间：

查看次数：

5455 次

最近记录：

5 年，1 月 前
相关归档
确定列中的NA值的数量 123
将MATLAB代码转换为R. 27
R中L的澄清 26
如何在R中停止/结束/暂停脚本？ 23
knitr显示整数的数字而没有科学记数法 15
Cowplot让ggplot2主题消失/如何查看当前的ggplot2主题,并恢复默认？ 15
Rmarkdown 中的 flextable 自动调整到 word doc 导致表格超出页边距 13
通过knitr,gWidgetsWWW使用时出现data.table错误 11
如何在累积总和图上的常规时间点获得插值的斜率？ 10
ggplot2 + gridExtra:如何确保不同大小的geog_bar绘图grob导致完全相同的条宽 10
难疑归档
如何更改远程Git存储库的URI(URL)？ 3529
如何在Git中获取当前分支名称？ 2321
finally块总是在Java中执行吗？ 2281
什么是NullReferenceException,我该如何解决？ 1876
使用jQuery将表单数据转换为JavaScript对象 1580
外部"C"在C++中有什么影响？ 1511
为什么要使用getter和setter/accessors？ 1472
如何从Git存储库中删除.DS_Store文件？ 1167
创建将T限制为枚举的通用方法 1122
Vim中的remap,noremap,nnoremap和vnoremap映射命令有什么区别？ 1045