R has some very useful libraries for working with spatial data. In this blog we will look at some of the libraries and demonstrate few basic functionalities. Lets start with reading a shapefile.
How to read a shapefile :
We will use the maptools package to read the shape file. Along with the maptools package, install the rgeos and sp packages. They will come handy later on.
To demonstrate reading a shapefile, we use the shapefile of US states which we download from here. The zip folder contains the file statesp020.shp which we will attempt to read. Lets first specify the projection that we want to use to read data in. The shapefile contains polygons in WGS84 projection, so lets define an object to hold that projection.
> library(maptools)
Loading required package: sp
Checking rgeos availability: TRUE
> crswgs84=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
The CRS class is present in the sp package and it holds the projection definition in proj4j format. We now read the shapefile
> states=readShapePoly("/home/mithil/java/R/statesp020.shp",proj4string=crswgs84,verbose=TRUE)
The shapefile is now read and stored into an object called “states”. readShapePoly function is used to read a shapefile that contains polygons. Lets explore the object that is created, beginning with what type it is.
> class(states) [1] "SpatialPolygonsDataFrame" attr(,"package") [1] "sp"
So the object is of type SpatialPolygonsDataFrame. This comes from the sp package. This object has 5 slots – data, polygons, bbox, plotOrder,bbox, proj4string. data contains the information about the polygons, polygons contains the actual polygon coordinates, bbox is the bounding box drawn around the boundaries of the shapefile and the proj4string is the projection.
> str(states@data)
'data.frame': 2895 obs. of 9 variables:
$ AREA : num 267 0 0 0 0 ...
$ PERIMETER : num 374.768 0.224 0.118 0.276 0.167 ...
$ STATESP020: int 2 3 4 5 6 7 8 9 10 11 ...
$ STATE : Factor w/ 53 levels "Alabama","Alaska",..: 2 2 2 2 2 2 2 2 2 2 ...
$ STATE_FIPS: Factor w/ 53 levels "01","02","04",..: 2 2 2 2 2 2 2 2 2 2 ...
$ ORDER_ADM : int 49 49 49 49 49 49 49 49 49 49 ...
$ MONTH_ADM : Factor w/ 12 levels "April","August",..: 5 5 5 5 5 5 5 5 5 5 ...
$ DAY_ADM : int 3 3 3 3 3 3 3 3 3 3 ...
$ YEAR_ADM : int 1959 1959 1959 1959 1959 1959 1959 1959 1959 1959 ...
- attr(*, "data_types")= chr "N" "N" "N" "C" ...
> str(states@polygons[[1]])
Formal class 'Polygons' [package "sp"] with 5 slots
..@ Polygons :List of 1
.. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. ..@ labpt : num [1:2] -152.7 64.6
.. .. .. ..@ area : num 267
.. .. .. ..@ hole : logi FALSE
.. .. .. ..@ ringDir: int 1
.. .. .. ..@ coords : num [1:70238, 1:2] -157 -157 -157 -157 -157 ...
..@ plotOrder: int 1
..@ labpt : num [1:2] -152.7 64.6
..@ ID : chr "0"
..@ area : num 267
> states@bbox
min max
x -179.13339 179.78821
y 17.67469 71.39805
> states@proj4string
CRS arguments:
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
How to plot a shapefile :
To see how the shapefile looks like or to create an image out of it use
> plot(states)
How to transform a shapefile :
To transform a shapefile to a different coordinate system use the spTransform method from the rgdal package. Lets transform our sp object to mercator projection
crsmerc=CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs")
states_transformed<-spTransform(states,CRS=crsmerc)
How to check if a point falls inside a polygon :
One last but useful function that we will see in this post is the gContains function from the rgeos package. This function checks whether a polygon contains a point or more generally whether a geometry contains another geometry.
> library(rgeos) > p<-SpatialPoints(list(lng,lat), proj4string=crswgs84) > gContains(fdg,pt)
gContains returns true or false depending on whether the polygon does or does not contain the point.
This post demonstrates the power of R in handling spatial data. The packages introduced in this post- maptools and rgeos are quite powerful and interested readers may want to go through the library documentation to see all the functionalities it provides