Building spatial objects with R

I started to work with geolocalized data some time ago. I had already performed basic spatial analysis but it was my (real) first steps in the wonderful world of Geographical Information Systems (GIS). I first used QGIS which is a pretty great (free!) tool to work with spatial data. However, I needed a more flexible tool and I decided to start learning how to perform spatial analysis with R. In this post I’m going to show you how to build and export spatial objects with R using the sp package. This package provides classes and methods for dealing with spatial data such as points, lines or polygons. We will focus on the creation of a spatial polygon that we will then export in shapefile format.

Let’s consider three 2D points p11, p12 and p13:

p11=c(2.1,40.1)
p12=c(2.1,40.2)
p13=c(2.2,40.15)

At the beginning, a polygon is nothing more than a matrix x1 containing (in row) the coordinates of the vertices of the polygon:

p1=rbind(p11,p12,p13)

This polygon can be plotted using the function polygon. The function polygon, like the function points, segments and lines, draws a shape on an existing plot. This is why we need to create an empty plot before drawing the polygon:

par(mar=c(0,0,0,0))
plot(1,type="n",axes=FALSE,xlab="",ylab="",xlim=c(min(p1[,1]),max(p1[,1])),ylim=c(min(p1[,2]),max(p1[,2])))
polygon(p1,col="steelblue3",border="#CC6666",lwd=5)

Then we transform our matrix p1 into a Polygon-class object P1 using the function Polygon:

P1=Polygon(rbind(p1,p1[1,]))

This new object P1 is a S4 object from the Spatial class, in which we can identify different components, called slots. To get the slot names of an S4 object we use the slotNames function and we need to use a @ to access a slot.

We can identify 5 slots:

  • labpt the centroid of the polygon,
  • area the surface of P1,
  • hole TRUE if the polygon is a hole,
  • ringDir the ring direction of the ring (polygon) coordinates, holes are expected to be anti-clockwise,
  • coords the coordinates of the polygon.

It is important to note that the first point in coords should equal the last point. Note that, as Thomas Louail (yes the famous LouBar!) has rightly pointed out to me, to transform the matrix p1 into a polygon, Polygon(p1) is enough, we do not need to “loop” the polygon manually since sp does it automatically.

The attribute hole allows us to create a hole inside a polygon, for example we can create a new polygon H1 which is (spatially) contained in P1 and defined as a hole:

h11=c(2.115,40.12)
h12=c(2.115,40.18)
h13=c(2.18,40.15)
h1=rbind(h11,h12,h13)
H1=Polygon(rbind(h1,h1[1,]),hole=TRUE)

We obtain this shape:

par(mar=c(0,0,0,0))
plot(1,type="n",axes=FALSE,xlab="",ylab="",xlim=c(min(p1[,1]),max(p1[,1])),ylim=c(min(p1[,2]),max(p1[,2])))
polygon(p1,col="steelblue3",border="#CC6666",lwd=5)
polygon(h1,col="green3",border="#CC6666",lwd=5)

Note that P1 and H1 are two distinct spatial objects. To define H1 as a hole of P1 we can merge them into one object using the function Polygons and create an object Ps1 which is a list of Polygon-class objects:

Ps1 = Polygons(list(P1,H1), ID = "P1")

This new object contains 5 slots:

  • the Polygons P1 and H1,
  • plotOrder order in which the Polygons should be plotted,
  • labpt the centroid,
  • the Polygons’ ID,
  • area the surface.

It seems that the centroid and the area provided in these slots are the area and the centroid of the first element P1 and not the ones of Ps1. To extract the area of H1, one can use the code below:

Ps1@Polygons[[2]]@area

More systematically, to obtain the area of each element of a Polygons we can use the function sapply:

sapply(Ps1@Polygons, function(x) x@area)

In the same way we can create another Polygons Ps2

p21=c(2.3,40.1)
p22=c(2.3,40.2)
p23=c(2.2,40.15)

p2=rbind(p21,p22,p23)
P2=Polygon(p2)
Ps2 = Polygons(list(P2), ID = "P2")

and merge Ps1 and Ps2 into a SpatialPolygons SPs:

SPs = SpatialPolygons(list(Ps1,Ps2),proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))

The main difference between a Polygons and a SpatialPolygons is the argument proj4string. This argument is a CRS class object that allows us to define the reference system of our object. In this case I used longitude/latitude coordinates. The reference is really important to manipulate spatial objects. This new object contains 4 slots:

  • the polygons Ps1 and Ps2,
  • plotOrder order in which the Polygons should be plotted,
  • bbox a bounding box that contains a matrix with the coordinate values,
  • proj4string a CRS class object, the coordinate system and projection information.

Finally we plot the results:

par(mar=c(0,0,0,0))
plot(SPs,col="steelblue3",border="#CC6666",lwd=5)

To export SPs in a shapefile format we can use the function writeOGR (package rgdal). But first we need to transform our SpatialPolygons into a SpatialPolygonsDataFrame:

data=data.frame(ID=c(1,2),Name=c("P1","P2"),row.names=c("P1","P2"))
SPs=SpatialPolygonsDataFrame(SPs,data)

A SpatialPolygonsDataFrame is a SpatialPolygons object with data attached (the attribute table) accessible through the slot data. One important rule, the row.names of the dataframe must match the SpatialPolygons’ IDs. We can finally export SPs using the function writeOGR.

writeOGR(obj=SPs,dsn="SPs",layer="SPs",driver="ESRI Shapefile")

The output is available in a folder SPs (dsn) located in your working directory. This folder is composed of four files:

  • SPs.shp the Polygons,
  • SPs.shx a positional index of the feature geometry,
  • SPs.dbf the data,
  • SPs.prj the coordinate system and projection information proj4string.

Finally, we can open the file with QGIS and check that everything is ok.

comments powered by Disqus