如何从 R 包“spatstat”转换样本数据集到形状文件中

发布于 2024-11-03 14:47:15 字数 550 浏览 2 评论 0原文

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

所以我尝试使用shapefiles 包的函数 convert.to.shapefile(4) 将 R 中 spatstat 包中包含的 bei 数据集转换为 shapefile。不幸的是,事实证明这比我想象的要困难。有人有这方面的经验吗?如果您愿意帮我一把,我将不胜感激。

谢谢, 瑞安

参考资料: spatstat, shapefile

I have written a kernel density estimator in Java that takes input in the form of ESRI shapefiles and outputs a GeoTIFF image of the estimated surface. To test this module I need an example shapefile, and for whatever reason I have been told to retrieve one from the sample data included in R. Problem is that none of the sample data is a shapefile...

So I'm trying to use the shapefiles package's funciton convert.to.shapefile(4) to convert the bei dataset included in the spatstat package in R to a shapefile. Unfortunately this is proving to be harder than I thought. Does anyone have any experience in doing this? If you'd be so kind as to lend me a hand here I'd greatly appreciate it.

Thanks,
Ryan

References:
spatstat,
shapefiles

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。

评论(2

终弃我 2024-11-10 14:47:15

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

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 的另一个参数,或者在使用 proj4string(spDF) <- CRS("+proj=etc. ..")(但这本身就是一个完整的问题,我们可以在其上编写页面)。

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

library(rgdal)

(在中使用writePolyShape如果rgdal不可用,则使用maptools包。

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

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

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

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

list.files(pattern = "^bei")

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

请注意,ppp 对象没有通用的“as.Spatial”转换器,因为决策必须确定这是否是一个要点带有标记的图案等等 - 尝试编写一个报告是否需要虚拟数据等可能会很有趣,

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

library(sp);小插图(“sp”)
库(spatstat);小插图(“spatstat”)

There are converter functions for Spatial objects in the spatstat and maptools packages that can be used for this. A shapefile consists of at least points (or lines or polygons) and attributes for each object.

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

Coerce bei to a Spatial object, here just points without attributes since there are no "marks" on the ppp object.

spPoints <- as(bei, "SpatialPoints")

A shapefile requires at least one column of attribute data, so create a dummy.

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

Using the SpatialPoints object and the dummy data, generate a SpatialPointsDataFrame.

spDF <- SpatialPointsDataFrame(spPoints, dummyData)

At this point you should definitely consider what the coordinate system used by bei is and whether you can represent that with a WKT CRS (well-known text coordinate reference system). You can assign that to the Spatial object as another argument to SpatialPointsDataFrame, or after create with proj4string(spDF) <- CRS("+proj=etc...") (but this is an entire problem all on its own that we could write pages on).

Load the rgdal package (this is the most general option as it supports many formats and uses the GDAL library, but may not be available because of system dependencies.

library(rgdal)

(Use writePolyShape in the maptools package if rgdal is not available).

The syntax is the object, then the "data source name" (here the current directory, this can be a full path to a .shp, or a folder), then the layer (for shapefiles the file name without the extension), and then the name of the output driver.

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

Note that the write would fail if the "bei.shp" already existed and so would have to be deleted first unlink("bei.shp").

List any files that start with "bei":

list.files(pattern = "^bei")

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

Note that there is no general "as.Spatial" converter for ppp objects, since decisions must be made as to whether this is a point patter with marks and so on - it might be interesting to try writing one, that reports on whether dummy data was required and so on.

See the following vignettes for further information and details on the differences between these data representations:

library(sp); vignette("sp")
library(spatstat); vignette("spatstat")

楠木可依 2024-11-10 14:47:15

一般的解决方案是:

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

例如,如果我们考虑 spatstat 中的 hamster 数据集:

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")

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

我没有对 bei 数据集执行上述操作的原因之一是它不包含任何数据,因此我们无法将其强制转换为 SpatialPointsDataFrame 对象。 我们可以使用的数据,在 bei.extra 中(与 bei 同时加载),但是这些额外的数据还是在常规网格上。因此,我们必须

  • bei.extra 转换为 SpatialGridDataFrame 对象(例如 bei.spg
  • SpatialPoints 对象(例如 bei.sp
  • overlay()bei.sp 指向 <代码>bei.spg grid,从网格中为 bei 中的每个点生成值,
  • 这应该为我们提供一个可以使用 writeOGR() 写出的 SpatialPointsDataFrame如上

如您所见,只是为了给您一个 Shapefile,这需要更多的工作。我展示的仓鼠数据示例是否足够?如果没有,我可以找到我的 Bivand 等明天并执行 bei 的步骤。

A general solution is:

  • convert the "ppp" or "owin" classed objects to appropriate classed objects from the sp package
  • use the writeOGR() function from package rgdal to write the Shapefile out

For example, if we consider the hamster data set from spatstat:

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

first convert this object to a SpatialPointsDataFrame object:

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

This gives us a sp object to work from:

> 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

This object has a single variable in the @data slot:

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

So say we now want to write out this variable as an ESRI Shapefile, we use writeOGR()

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

This will create several marks.xxx files in directory hamster created in the current working directory. That set of files is the ShapeFile.

One of the reasons why I didn't do the above with the bei data set is that it doesn't contain any data and thus we can't coerce it to a SpatialPointsDataFrame object. There are data we could use, in bei.extra (loaded at same time as bei), but these extra data or on a regular grid. So we'd have to

  • convert bei.extra to a SpatialGridDataFrame object (say bei.spg)
  • convert bei to a SpatialPoints object (say bei.sp)
  • overlay() the bei.sp points on to the bei.spg grid, yielding values from the grid for each of the points in bei
  • that should give us a SpatialPointsDataFrame that can be written out using writeOGR() as above

As you see, that is a bit more involved just to give you a Shapefile. Will the hamster data example I show suffice? If not, I can hunt out my Bivand et al tomorrow and run through the steps for bei.

~没有更多了~
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文