How to convert a sample dataset from the R package "spatstat" into a shapefile

Question

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

Answer 1

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

Answer 2

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.