在R中创建Shapefile

Question

我正在尝试在R中创建一个形状文件，之后将其导入到Fusion Table或其他GIS应用程序中。

首先，我导入了一个包含加拿大所有普查区的空白形状文件。 我根据CT的唯一ID附加了其他数据（以表格格式），并绘制了我的结果。 目前，我只需要温哥华的数据，并希望导出一个仅包含温哥华CT和我新附加的属性数据的形状文件。

以下是我的代码（由于隐私原因省略了某些部分）：

shape <- readShapePoly('C:/TEST/blank_ct.shp') #Load blank shapefile
shape@data = data.frame(shape@data, data2[match(shape@data$CTUID, data2$CTUID),]) #data2 is my created attributes that I'm attaching to blank file
shape1 <-shape[shape$CMAUID == 933,] #selecting the Vancouver CTs

我看到其他人使用这个: writePolyShape 来创建shapefile。我尝试过，并且它在某种程度上是有效的。它创建了 .shp，.dbf 和 .shx 文件。但我缺少 .prj 文件，也不确定如何创建它。是否有更好的方法来创建shapefile呢？

如果您能提供任何帮助，将不胜感激。

- user2142810

---

这是我写文件的方式：writeOGR(obj = opno.skupaj.mean[[1]], dsn = "q:/path/to/file/spat_odstrel_skupaj.shp", layer = "spat_odstrel_skupaj", driver = "ESRI Shapefile")。请注意，图层和文件名是相同的（去掉了 .shp）。 - Roman Luštrik

Answer 1

使用rgdal和writeOGR。rgdal将保留投影信息。

类似这样的内容。

library(rdgal)

shape <- readOGR(dsn = 'C:/TEST', layer = 'blank_ct')
# do your processing
shape@data = data.frame(shape@data, data2[match(shape@data$CTUID, data2$CTUID),]) #data2 is my      created attributes that I'm attaching to blank file
shape1 <-shape[shape$CMAUID == 933,]
writeOGR(shape1, dsn = 'C:/TEST', layer ='newstuff', driver = 'ESRI Shapefile')

请注意，dsn是包含.shp文件的文件夹，而layer是不带.shp扩展名的shapefile名称。它将读取（readOGR）和写入（writeOGR）所有组件文件（.dbf，.shp，.prj等）。

- mnel

---

这样做不行，你没有在writeOGR调用中指定shape或shape1对象 - 它怎么知道要写什么？ - Spacedman

我认为 writeOGR 的 dsn 路径需要是一个以 .shp 结尾的完整文件路径。 - Roman Luštrik

@Roman 不需要那个，可以使用writeOGR(x, ".", "layername", "ESRI Shapefile")或writeOGR(x, "C:/temp/layername.shp", "layername", "ESRI Shapefile")。 - mdsumner

@mdsumner 哈！我想知道我做错了什么。我以前用你的方法写不出来，但现在可以了。我们是不是跨越了奇怪事情经常发生的区域…… - Roman Luštrik

@mdsumner很重要的一点是，如果决定不指定完整文件路径，则dsn中的路径不包括尾部斜杠。 - Roman Luštrik

Answer 2

问题解决了！再次感谢那些帮助我的人！

这是我最终所做的：

按照 Mnel 所写的，这条命令将创建一个 shapefile。

writeOGR(shape1, dsn = 'C:/TEST', layer ='newstuff', driver = 'ESRI Shapefile')

然而，当我运行这行代码时，出现了以下错误：

Can't convert columns of class: AsIs; column names: ct2,mprop,mlot,mliv

这是因为我的属性数据不是数字，而是字符。幸运的是，我的属性数据都是数字，所以我运行了transform()来解决这个问题。

shape2 <-shape1
shape2@data <- transform(shape1@data, ct2 = as.numeric(ct2),
mprop = as.numeric(mprop),
mlot = as.numeric(mlot),
mliv = as.numeric(mliv))

我再次尝试使用writeOGR()命令，但我仍未获得我要找的.prj文件。问题在于我在导入文件时没有指定shapefile的坐标系。既然我已经知道坐标系是什么，那么我所要做的就是在导入时定义它。

readShapePoly('C:/TEST/blank_ct.shp',proj4string=CRS("+proj=longlat +datum=WGS84")

之后，我重新运行了我想要使用shapefile进行的所有操作，包括用于导出的writeOGR命令。就是这样了！

- user2142810