Различия проекции в R с использованием sf и sp

У меня есть сетка, которую я преобразовал из GeoTIFF в шейп-файл. Я хотел бы преобразовать и экспортировать шейп-файл как GeoPackage и изменить проекцию, чтобы она использовала Британскую национальную сетку в качестве географической системы координат при открытии в ГИС. Однако, похоже, это работает только с 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 sf sp projection coordinate-systems
person Chris    schedule 27.07.2018    source источник

Ответы (3)


arrow_upward
1
arrow_downward

Решение этой проблемы (благодаря тому, что Роджер опубликовал здесь), использует lwgeom, чтобы выполнить преобразование. В сообщении Роджера на 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")
person Chris    schedule 03.08.2018

arrow_upward
0
arrow_downward

Используя 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?

person Edzer Pebesma    schedule 31.07.2018
comment
Спасибо, что посмотрели на Эдзер. Однако я не думаю, что это проблема. В примерах v1 и v2 я тестировал различные способы определения crs с помощью st_transform. v1 использует код epsg, а v2 - это то место, где я использовал проекцию, которая, как я знаю, содержит данные. Однако датум не отображается, когда вы пытаетесь использовать st_crs (BNG_Grid_sf_v1) или st_crs (BNG_Grid_sf_v2). Сравнивая это с sp-эквивалентом proj4string (BNG_Grid_sp), у которого есть датум, можно предположить, что проблема может заключаться в том, что датум для BNG не сохраняется sf при использовании st_transform, и ArcMap не может использовать то, чего там нет. - person Chris; 01.08.2018

arrow_upward
0
arrow_downward

Следуйте за https://github.com/r-spatial/sf/issues/810, здесь нет.

person Roger Bivand    schedule 01.08.2018