为什么 sf::st_transform() 返回的对象与调用中使用的投影不同？
Why is sf::st_transform() returning an object with a different projection than used in the call?
我想使用 sf::st_transform() 重新投影我的 sf 对象，但转换对象的投影与我在转换调用中指定的投影不同。为什么？
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"
投影之间的唯一区别是 proj4string 中的 +x_0 元素。 p2 的投影中删除了尾随的 1e-10。
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
我对制图学的了解不够，无法判断差异是由于帮助页面中提到的逆投影问题还是其他原因造成的。
相关：
我想使用 sf::st_transform() 重新投影我的 sf 对象，但转换对象的投影与我在转换调用中指定的投影不同。为什么？
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"
投影之间的唯一区别是 proj4string 中的 +x_0 元素。 p2 的投影中删除了尾随的 1e-10。
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
我对制图学的了解不够，无法判断差异是由于帮助页面中提到的逆投影问题还是其他原因造成的。
相关：