Почему sf :: st_transform () возвращает объект с проекцией, отличной от используемой в вызове? - PullRequest
Новая
       98

Почему sf :: st_transform () возвращает объект с проекцией, отличной от используемой в вызове?

0 голосов
/ 06 апреля 2019

Я хочу перепроектировать мой sf-объект, используя sf::st_transform(), но проекция преобразованного объекта отличается от проекции, которую я указал в вызове преобразования. Почему? 

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

# target proj4string
my_crs <- "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"

# create test sfc
p1 <- st_sfc(st_point(c(1,2)), crs = "+init=epsg:3857")

# re-project p1 to `my_crs`
p2 <- st_transform(p1, crs = my_crs)

all.equal(my_crs, st_crs(p2)$proj4string)
#> [1] "1 string mismatch"

st_crs(p2)
#> Coordinate Reference System:
#>   No EPSG code
#>   proj4string: "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"

Единственная разница между проекциями - это элемент +x_0 в строке proj4. Трейлинг 1e-10 был удален в проекции p2.

1 Ответ

0 голосов
/ 06 апреля 2019

sf::st_transform() использует GDAL для проекций. Страница справки гласит: 

Some PROJ.4 projections are not supported by GDAL, e.g. "+proj=wintri" because it does 
not have an inverse projection. Projecting to unsupported projections can be done by 
st_transform_proj, part of package lwgeom. Note that the unsupported proj4string cannot 
be passed as argument to st_crs, but has to be given as character string.

Если вы используете PROJ.4 для проекции в вашем примере, она работает правильно: 

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(lwgeom)
#> Linking to liblwgeom 2.5.0dev r16016, GEOS 3.6.1, PROJ 4.9.3

my_crs <- "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"

p1 <- st_sfc(st_point(c(1,2)), crs = "+init=epsg:3857")

p2 <- lwgeom::st_transform_proj(p1, crs = my_crs)

all.equal(my_crs, st_crs(p2)$proj4string)
#> [1] TRUE

Я недостаточно знаю о картографии, чтобы сказать, связано ли это с проблемой обратной проекции, упомянутой на странице справки, или с чем-то еще.

Связанный:

https://stackoverflow.com/a/51663647/1707525

https://github.com/r-spatial/sf/issues/810

https://github.com/r-spatial/sf/issues/1019

...