6

我用 Java 编写了一个内核密度估计器，它以 ESRI shapefile 的形式输入并输出估计表面的 GeoTIFF 图像。为了测试这个模块，我需要一个示例 shapefile，无论出于何种原因，我都被告知要从 R 中包含的示例数据中检索一个。问题是所有示例数据都不是 shapefile...

所以我正在尝试使用 shapefiles 包的convert.to.shapefile(4)功能将 R 中 spatstat 包中包含的 bei 数据集转换为 shapefile。不幸的是，事实证明这比我想象的要难。有没有人有这方面的经验？如果你能在这儿帮我一把，我将不胜感激。

谢谢，瑞安

参考： spatstat , shapefiles

4

2 回答 2

6

和包中的Spatial对象有转换器函数可用于此目的。shapefile 至少由每个对象的点（或线或多边形）和属性组成。spatstatmaptools

library(spatstat)
library(sp)
library(maptools)
data(bei)

强制bei到一个Spatial对象，这里只是没有属性的点，因为对象上没有“标记” ppp。

spPoints <- as(bei, "SpatialPoints")

shapefile 至少需要一列属性数据，因此请创建一个虚拟文件。

dummyData <- data.frame(dummy = rep(0, npoints(bei)))

使用SpatialPoints对象和虚拟数据，生成一个SpatialPointsDataFrame.

spDF <- SpatialPointsDataFrame(spPoints, dummyData)

在这一点上，您绝对应该考虑使用的坐标系bei是什么以及您是否可以用WKT CRS（众所周知的文本坐标参考系）来表示它。Spatial您可以将它作为另一个参数分配给对象SpatialPointsDataFrame，或者在 create with 之后分配proj4string(spDF) <- CRS("+proj=etc...")（但这是一个完全独立的问题，我们可以在上面写页面）。

加载rgdal包（这是最通用的选项，因为它支持多种格式并使用 GDAL 库，但由于系统依赖性可能不可用。

library(rgdal)

（如果没有，请writePolyShape在maptools包装中使用）。rgdal

语法是对象，然后是“数据源名称”（这里是当前目录，可以是 .shp 或文件夹的完整路径），然后是层（对于 shapefile，是不带扩展名的文件名），以及然后是输出驱动程序的名称。

writeOGR(obj = spDF, dsn = ".", layer = "bei", driver = "ESRI Shapefile")

请注意，如果“bei.shp”已经存在，则写入将失败，因此必须先删除unlink("bei.shp")。

列出所有以“bei”开头的文件：

list.files(pattern = "^bei")

[1] "bei.dbf" "bei.shp" "bei.shx"

请注意，没有用于ppp对象的通用“as.Spatial”转换器，因为必须决定这是否是带有标记的点模式等等 - 尝试编写一个报告虚拟数据是否为需要等等。

有关这些数据表示之间差异的更多信息和详细信息，请参阅以下小插图：

图书馆（SP）；小插图（“sp”）库（spatstat）；小插图（“spatstat”）

于 2011-04-26T21:32:52.220 回答
3

一个通用的解决方案是：

  • 将包中的"ppp"或"owin"分类对象转换为适当的分类对象sp
  • 使用writeOGR()包中的函数rgdal写出Shapefile

例如，如果我们考虑hamster来自 的数据集spatstat：

require(spatstat)
require(maptools)
require(sp)
require(rgdal)
data(hamster)

首先将此对象转换为SpatialPointsDataFrame对象：

ham.sp <- as.SpatialPointsDataFrame.ppp(hamster)

这给了我们一个sp工作对象：

> str(ham.sp, max = 2)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 303 obs. of  1 variable:
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:303, 1:2] 6 10.8 25.8 26.8 32.5 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..@ bbox       : num [1:2, 1:2] 0 0 250 250
  .. ..- attr(*, "dimnames")=List of 2
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots

@data该对象在槽中有一个变量：

> head(ham.sp@data)
     marks
1 dividing
2 dividing
3 dividing
4 dividing
5 dividing
6 dividing

所以说我们现在想把这个变量写成一个 ESRI Shapefile，我们使用writeOGR()

writeOGR(ham.sp, "hamster", "marks", driver = "ESRI Shapefile")

这将在当前工作目录中创建的目录中创建几个marks.xxx文件。hamster那组文件就是ShapeFile。

我没有对bei数据集执行上述操作的原因之一是它不包含任何数据，因此我们无法将其强制为SpatialPointsDataFrame对象。我们可以在（与 加载的同时）中使用一些数据，但这些额外的数据或在常规网格上。所以我们不得不bei.extrabei

  • 转换bei.extra为SpatialGridDataFrame对象（例如bei.spg）
  • 转换bei为SpatialPoints对象（例如bei.sp）
  • overlay()网格上的bei.sp点bei.spg，从网格中为每个点生成值bei
  • 那应该给我们一个SpatialPointsDataFrame可以用writeOGR()上面写出来的

如您所见，这只是为了给您一个 Shapefile。我展示的hamster数据示例是否足够？如果没有，我明天可以找到我的Bivand 等人bei，然后完成.

于 2011-04-26T21:28:32.037 回答