使用 sf 和 sp 在 R 中的投影差异答案

【问题标题】：Projection differences in R using sf and sp使用 sf 和 sp 在 R 中的投影差异
【发布时间】：2018-07-27 08:24:06
【问题描述】：

我有一个从 GeoTIFF 转换为 shapefile 的网格。我想将 shapefile 转换并导出为 GeoPackage 并更改投影，以便在 GIS 中打开时使用英国国家网格作为地理坐标系。然而，这似乎只适用于使用sp 而不是sf（它似乎没有保留基准等方面）。

这是一个问题，因为我想导出包含多个图层的 GeoPackages，您目前只能在 sf 而不是 sp 中执行此操作。我做错了吗？

library(rgdal)
library(sf)

download.file("https://drive.google.com/uc?id=1URbux7Sw25KFTySqRFKXk53DV2UK4lsA&export=download" , destfile="Grid_Shapefile.zip")
unzip("Grid_Shapefile.zip")
Grid_sp <- readOGR(".", "Grid_Shapefile")
Grid_sf <- st_as_sf(Grid_sp)

BNG_Grid_sp <- spTransform(Grid_sp, CRS("+init=epsg:27700"))
BNG_Grid_sf_v1 <- st_transform(Grid_sf, crs=27700)
BNG_Grid_sf_v2 <- st_transform(Grid_sf, crs="+init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894")

BNG_Grid_sf_v1_geom <- st_geometry(BNG_Grid_sf_v1)
BNG_Grid_sf_v2_geom <- st_geometry(BNG_Grid_sf_v2)

proj4string(BNG_Grid_sp)
attributes(BNG_Grid_sf_v1_geom)
attributes(BNG_Grid_sf_v2_geom)

writeOGR(BNG_Grid_sp, dsn = "BNG_Grid_sp.gpkg", layer = "Grid_sp", driver = "GPKG")
st_write(BNG_Grid_sf_v1, "BNG_Grid_sf_v1.gpkg", "Grid_sf_v1")
st_write(BNG_Grid_sf_v2, "BNG_Grid_sf_v2.gpkg", "Grid_sf_v2")

【问题讨论】：

标签： r projection coordinate-systems sp sf

【解决方案1】：

解决方案（感谢 Roger 发布的 here）是使用 lwgeom 包进行转换。Roger 在 sf GitHub 上的帖子提供了更多详细信息。

library(rgdal)
library(sf)

download.file("https://drive.google.com/uc?id=1URbux7Sw25KFTySqRFKXk53DV2UK4lsA&export=download" , destfile="Grid_Shapefile.zip")
unzip("Grid_Shapefile.zip")
Grid_sp <- readOGR(".", "Grid_Shapefile")
Grid_sf <- st_as_sf(Grid_sp)

library(lwgeom)
BNG_Grid_sf_v4 <- st_transform_proj(Grid_sf, crs="+init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894")
st_crs(BNG_Grid_sf_v4)
st_write(BNG_Grid_sf_v4, "BNG_Grid_sf_v4.gpkg", "Grid_sf_v4")

【讨论】：

【解决方案2】：

使用ogrinfo，尤其是命令

ogrinfo BNG*v1.gpkg Grid_sf_v1 > info1
ogrinfo BNG*v2.gpkg Grid_sf_v1 > info2

两个 info[1|2] 文件之间的差异，除了明显的 _v1 _v2 命名之外，还有：

13c13
<             TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489]],
---
>             TOWGS84[446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894]],

_v2 中的额外数字是否会导致您在 ArcGIS 中遇到麻烦？

【讨论】：

感谢您查看此 Edzer。但是，我认为这不是问题所在。 v1 和 v2 示例是我测试使用 st_transform 定义 crs 的不同方法。 v1 使用 epsg 代码，而 v2 是我使用我知道包含基准的投影的地方。但是，当您尝试 st_crs(BNG_Grid_sf_v1) 或 st_crs(BNG_Grid_sf_v2) 时，基准不会出现。将此与确实具有基准的 sp 等效 proj4string(BNG_Grid_sp) 进行比较，表明问题可能是使用 st_transform 时 sf 未保留 BNG 的基准，并且 ArcMap 无法使用不存在的东西。跨度>

【解决方案3】：

请关注https://github.com/r-spatial/sf/issues/810，而不是这里。

【讨论】：