将缺失的CRS添加到R简单特征中的sf对象

问题描述 投票:0回答:2

有没有办法添加缺失的CRS?我收到了一些旧的 Esri shapefile,其中很大一部分有

CRS: NA
,我希望能够批量设置丢失的 CRS。当我导入QGIS时我可以手动设置它,但似乎无法使用R来设置它。我尝试了
st_transform()
,但如果一开始就没有设置CRS,它似乎不起作用。

r mapping gis r-sf
2个回答
0
投票

感谢您的评论。我尝试了

st_crs <-
st_set_crs()
但似乎仍然无法正确设置或绘制形状文件。我注意到那些缺少
.prj
文件的文件是没有 CRS 的文件,这是有道理的。所以我写了一个块来为目录中的每个 shapefile 创建一个 proj 文件(如果不存在的话)。有点破解,但它有效:)

编辑:添加了一个粗略的日志系统来跟踪更改

p = "/path/to/shape/files"

proj_str_bng = 'PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.999601272],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0]]'

proj_str_wgs84 = 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'

shp_to_prj <- list.files(paste0(wd, p), pattern = "\\.shp$")

if (!file.exists(paste0(wd,p,"prj_log"))){
write("file, dt", file = paste0(wd,p,"prj_log"))
}

for (i in 1:length(shp_to_prj)) {
  shp_to_prj[[i]] <- paste0(wd, p, str_remove_all(shp_to_prj[[i]], ".shp"), ".prj")
  if(file.exists(shp_to_prj[[i]])) {
    next
  }
  write(proj_str_bng, file = shp_to_prj[[i]])
  write(paste0( shp_to_prj[[i]], ",",Sys.time()), file = paste0(wd,p,"prj_log"), append=TRUE)
}
write(paste0( "^===============================================================^", ",",Sys.time()), file = paste0(wd,p,"prj_log"), append=TRUE)


0
投票

我知道您提到尝试使用 st_set_crs() 但您需要确保您使用的是非常具体的语法(参见here)我试图使用普通语法来使用 st_set_crs(),例如。

st_set_crs(myObj,crs)

直到使用“管道语法”才起作用

myObj = myObj %>% st_set_crs(myEPSGNumber)
.

我认为你也可以在最初读取shapefile时设置crs:

myObj = st_read("myObj.shp",crs=st_crs(myEPSGNumber)

我花了一段时间才弄清楚,但希望能帮助你或其他任性的灵魂!

© www.soinside.com 2019 - 2024. All rights reserved.