R and GIS – working with shapefiles

April 18, 2016
By

(This article was first published on R – StudyTrails, and kindly contributed to R-bloggers)

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([email protected])
'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([email protected][[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

> [email protected]
         min       max
x -179.13339 179.78821
y   17.67469  71.39805
> [email protected]
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)

Selection_172

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 [email protected] +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

To leave a comment for the author, please follow the link and comment on their blog: R – StudyTrails.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)