无法投影简单特征或更改投影

问题描述

我正在尝试将csv转换为sf空间数据文件，但是我遇到了我无法弄清的错误。

示例：

library(tidyverse)
library(sf)
#> Linking to GEOS 3.8.0,GDAL 3.0.4,PROJ 6.3.1

point_df <- tibble::tribble(
     ~city_name,~longitude,~latitude,"Akron",-81.5190053,41.0814447,"Albany",-73.7562317,42.6525793,"Schenectady",-73.9395687,42.8142432,"Albuquerque",-106.650422,35.0843859,"Allentown",-75.4714098,40.6022939,"Bethlehem",-75.3704579,40.6259316,"Atlanta",-84.3879824,33.7489954,"Augusta",-82.0105148,33.4734978,"Austin",-97.7430608,30.267153,"Bakersfield",-119.0187125,35.3732921
  )

point_sf <- st_as_sf(point_df,coords = c("longitude","latitude"))

point_sf <- st_set_crs(point_sf,4326)

st_transform(point_sf,102003)
#> Warning in CPL_crs_from_input(x): GDAL Error 1: PROJ: proj_create_from_database:
#> crs not found
#> Error in CPL_transform(x,crs,aoi,pipeline,reverse): crs not found: is it missing?

任何帮助将不胜感激。

编辑我找到了一个适用于github页面的kludgy解决方案，但如果可能的话，我仍在寻找更系统的解决方案。 https://github.com/r-spatial/sf/issues/1419

这里的解决方案是将sf对象转换为sp，然后再更改回sf。

reProject <- function (sf,proj_in = "+init=epsg:4326",proj_out = "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs") {
  require(sp)
  
  data_sp <- as(sf,"Spatial")
  
  proj4string(data_sp) <- CRS(proj_in)
  
  sf_out <- st_as_sf(spTransform(data_sp,CRS(proj_out)))
}

dat_out <- reProject(point_sf)

解决方法

以下代码行似乎发生了某些事情。但是那没有发生。

point_sf <- st_as_sf(point_df,coords = c("longitude","latitude"))

尽管此行代码创建简单要素几何点对象，但此代码不创建简单要素几何列（sfc）对象。并且由于没有sfc对象，因此下一行代码不起作用。

point_sf <- st_set_crs(point_sf,4326)

在另一行代码中，函数st_set_crs（）从sf或sfc对象检索坐标参考系统。但是sf或sfc对象目前都不存在。

因此，必须先创建sfc对象，然后才能使用功能：st_set_crs（）。

在进行这些类型的简单要素项目时，确实可以按照以下步骤进行操作。

x.sfg <- st_multipoint(c(lon,lat),dim = "XY")   # create sf geometry from lon/lat
x.sfc <- st_sfc(x.sfg,crs = 4326)                # create sfc from geometry 
x.sf <- st_sf(df,x.sfc)                          # create sf object from sfc

首先将对数和纬度转换为向量，然后创建矩阵，然后以正确的顺序创建简单要素对象。

lon <- c(-81.5190053,-73.7562317,-73.9395687,-106.650422,-75.4714098,-75.3704579,-84.3879824,-82.0105148,-97.7430608,-119.0187125) 
lat <- c(41.0814447,42.6525793,42.8142432,35.0843859,40.6022939,40.6259316,33.7489954,33.4734978,30.267153,35.3732921)
m <- matrix(data = c(lon,nrow = 10,ncol = 2,byrow  = FALSE)
m.sfg <- st_multipoint(m,dim = "XY")
m.sfc <- st_sfc(m.sfg,crs = 4326)
m.sf <- st_sf(df,m.sfc)
head(m.sf,3)

然后创建美国大陆的基本图，然后将简单要素对象绘制到基本图上。

plot(US_48,axes = TRUE)
plot(m.sf,add= TRUE,pch = 19,col = "red")

上面显示的与问题的链接似乎与此问题没有任何联系。此处显示的答案不会将sf对象转换为sp，然后又更改回sf。

该图显示在link：

r r sf spatial