1

所以我有两个 SF 对象，一个多边形和一个多边形类，代表一个城市的社区和子社区。两者都来自 ArcGIS，并且覆盖完全相同的区域。但是，当我绘制这两个对象时，子社区在我的传单图中略有错位， 如此处所示。原来的预测是：

iver `ESRI Shapefile'
Simple feature collection with 92 features and 19 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 55500 ymin: 428647.4 xmax: 101032.6 ymax: 447000
epsg (SRID):    NA
proj4string:    +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs

ESRI Shapefile'
Simple feature collection with 7 features and 7 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 438682.1 ymin: 6771629 xmax: 512270.9 ymax: 6800944
epsg (SRID):    3857
proj4string:    +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs

之后，我使用 st_transform 将两者都更改为 WGS84。

传单代码是：

   leaflet(width = "100%") %>%
  addProviderTiles("Stamen.Terrain") %>%
  setView(lng = 4.314960, lat = 51.916024, zoom = 10) %>%
  addPolygons(data = rtdm_gebieden, weight = 2, color = "black", fillOpacity = 0.8, fillColor = groen_kleuren) %>%
  addPolygons(data = rtdm_buurten, weight = 2, color = "red")

我怀疑这是因为 sf 对象的 bbox 略有不同

SF 对象 1

Geometry set for 92 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 3.940974 ymin: 51.84307 xmax: 4.602129 ymax: 52.0055
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

科幻对象 2

Geometry set for 7 features 
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 3.940748 ymin: 51.84212 xmax: 4.601808 ymax: 52.00453
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

我尝试将一个的 bbox 设置为另一个的 bbox，但使用 st_bbox 不起作用。任何帮助将不胜感激！

4

1 回答 1

1

在评论中讨论后找到的解决方案

不同空间数据集之间的这种轻微变化似乎很常见，尤其是在比利时（例如比利时兰伯特 EPSG 31370）和荷兰空间层中。

第一个 SF 对象在导入时显示：

epsg (SRID):    NA
proj4string:    +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs

EPSG 代码丢失，对 proj4 字符串开头的快速搜索指向 Amersfoort 投影 EPSG:28992，并带有以下 proj4 字符串：

+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs

+towgs84在原始 SF 对象中的末端部分是观察到的偏移的来源。取决于参考数据库，它是否存在......
根据这个 Postgis 论坛答案的 Jan Hartman 所说：

问题是国家格网的基准与 WGS84 不同：它使用稍微不同的椭球来近似地球表面。使用标准 PROJ epsg 文件，将不会计算此基准差异，并且您会得到最大 100 米的差异。要完全转换到 WGS84，您需要在 PROJ 参数字符串中添加一个额外的“+towgs”参数。

@Sjoerd Braaksma 应用的一种解决方案是返回生成图层的软件（即ArcGIS）并修复那里的投影。

另一种解决方案是，当您在 R 中导入数据时，强制 R 通过 epsg 代码分配投影，以确保使用正确的投影：

SF <- read_sf(your_layer, crs = 28992)

这将导致以下警告，此处可以忽略，因为我们确实不想重新投影此数据：

Warning message:
st_crs<- : replacing crs does not reproject data; use st_transform for that  warning

这是一个比利时数据集的示例：

> SF <- read_sf(dsn = "mylayer.shp")
> SF
## Simple feature collection with 2664 features and 19 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 22265.45 ymin: 21162.99 xmax: 295157.4 ymax: 244027.9
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49.8333339 +lat_2=51.16666733333333 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.01256 +y_0=5400088.4378 +ellps=intl +units=m +no_defs

强制投影：

> SF <- read_sf(dsn = "mylayer.shp", crs = 31370)
## Warning message:
## st_crs<- : replacing crs does not reproject data; use st_transform for that 

> SF
## Simple feature collection with 2664 features and 19 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 22265.45 ymin: 21162.99 xmax: 295157.4 ymax: 244027.9
## epsg (SRID):    31370
## proj4string:    +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs
于 2018-02-15T21:49:51.023 回答