Spatial Data

Vector Data

Lots of spatial data classes

  • From the sp package
  • Based on data frames
  • Each row has an associated geometry
  • Points, Lines, Polygons, Grids
  • Plot, reproject, overlay...
  • Read/Write based on GDAL/OGR library

Raster Data

From the raster package

  • A bit like matrices or arrays
  • Single or multiple layers
  • Numbers or categories
  • Analysis functions
  • Read/Write based on GDAL/OGR library

Example

I always find myself using these...

> require(sp)
> require(raster)
> require(rgdal)

Background Maps

The dismo package can get a Google Map tile

> require(dismo)
> map = gmap("Monterey Bay", exp = 2)
> plot(map)
plot of chunk unnamed-chunk-2
> class(map)
[1] "RasterLayer"
attr(,"package")
[1] "raster"

But google map data has restrictive re-use license.

OpenStreetMap

A crowd-sourced and open-licensed world map. With map tiles:

> require(OpenStreetMap)
Loading required package: OpenStreetMap
Loading required package: rJava
> monterey = openmap(c(37.135, -122.343), 
     c(36.433, -121.628))
> plot(monterey)
plot of chunk unnamed-chunk-3

Read a Shapefile

A common format for geographical data

  • Point, line, or polygon geometry for each feature
  • Attributes for each feature
  • Open standard - reads into any GIS
  • One shapefile is multiple actual files...
Data/monterey_otters.dbf
Data/monterey_otters.prj
Data/monterey_otters.shp
Data/monterey_otters.shx


Read in with readOGR

> mbcensus = readOGR("./Data", "monterey_otters")
OGR data source with driver: ESRI Shapefile 
Source: "./Data", layer: "monterey_otters"
with 569 features and 16 fields
Feature type: wkbPolygon with 2 dimensions
> class(mbcensus)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
> summary(mbcensus)
Object of class SpatialPolygonsDataFrame
Coordinates:
      min     max
x -208762 -154874
y -174277  -97353
Is projected: TRUE 
proj4string :
[+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120
+x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
+ellps=GRS80 +towgs84=0,0,0]
Data attributes:
   PERIMETER          AREA                 DEPTH    
 Min.   : 1.01   Min.   :0.056   -30 to -60m  :283  
 1st Qu.: 3.27   1st Qu.:0.431   0 to -30m    :285  
 Median : 4.76   Median :0.762   0 to -30m bay:  1  
 Mean   : 5.66   Mean   :0.947                      
 3rd Qu.: 6.89   3rd Qu.:1.276                      
 Max.   :39.07   Max.   :3.840                      
                                                    
    ATOS_ID        HAB_ID        ACRES      
 Min.   :180   180n   :  1   Min.   : 13.7  
 1st Qu.:252   181n   :  1   1st Qu.:106.6  
 Median :322   182n   :  1   Median :188.3  
 Mean   :322   182o   :  1   Mean   :234.0  
 3rd Qu.:393   183n   :  1   3rd Qu.:315.3  
 Max.   :464   183o   :  1   Max.   :948.9  
               (Other):563                  
    HECTARES     ZONE_CODE      ZONE          Year     
 Min.   :  5.6   b:  1     Min.   :1.0   Min.   :2012  
 1st Qu.: 43.1   n:285     1st Qu.:2.0   1st Qu.:2012  
 Median : 76.2   o:283     Median :2.0   Median :2012  
 Mean   : 94.7             Mean   :2.5   Mean   :2012  
 3rd Qu.:127.6             3rd Qu.:3.0   3rd Qu.:2012  
 Max.   :384.0             Max.   :3.0   Max.   :2012  
                                                       
    POLY_ID       dens_sm          pupratio     
 180n   :  1   Min.   : 0.000   Min.   :0.0000  
 181n   :  1   1st Qu.: 0.251   1st Qu.:0.0677  
 182n   :  1   Median : 0.654   Median :0.1367  
 182o   :  1   Mean   : 2.802   Mean   :0.1290  
 183n   :  1   3rd Qu.: 2.016   3rd Qu.:0.1905  
 183o   :  1   Max.   :19.857   Max.   :0.2979  
 (Other):563                                    
    lin_dens        trend5yr           Sect_ID   
 Min.   :0.111   Min.   :-0.05574   Min.   :2.0  
 1st Qu.:1.238   1st Qu.:-0.01155   1st Qu.:2.0  
 Median :3.048   Median : 0.00248   Median :3.0  
 Mean   :3.129   Mean   : 0.00667   Mean   :3.1  
 3rd Qu.:4.794   3rd Qu.: 0.01941   3rd Qu.:4.0  
 Max.   :6.968   Max.   : 0.09217   Max.   :5.0  
                                                 
> hist(mbcensus$pupratio)
plot of chunk unnamed-chunk-5
> plot(mbcensus)
plot of chunk unnamed-chunk-5

Coordinate Systems

Geographical data uses different coordinate systems

  • Lat-long
  • Ordnance Survey National Grid
  • Google Spherical Mercator
  • UTM zones
  • thousands more...

Luckily there's a standard way of specifying them, and a repository of short codes

  • +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs
  • +init=epsg:27700
  • +init=epsg:4326
  • +init=epsg:3857

Convert between them with spTransform

Transform

Get coordinate systems, transform

> proj4string(mbcensus)
[1] "+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
> proj4string(map)
[1] "+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"
> osm()
CRS arguments:
 +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 
> mbcensus = spTransform(mbcensus, 
     osm())

Mapping

Now I can overlay census regions on the map

> plot(monterey)
> plot(mbcensus, lty = 3, add = TRUE)
plot of chunk unnamed-chunk-7

Colour by trend

> require(classInt)
> ci = classIntervals(mbcensus$trend5yr, 
     10, style = "pretty")
> plot(monterey)
> colours = findColours(ci, c("red", 
     "grey", "blue"))
> plot(mbcensus, col = colours, border = NA, 
     add = TRUE)
> legend("bottomleft", title = "5-year trend", 
     fill = attr(colours, "palette"), 
     legend = names(attr(colours, 
         "table")), bg = "white")
plot of chunk unnamed-chunk-8

Point Data

Hydrological survey data from NOAA

survey_id,lat,long,depth,quality_code,active
H05247,36.895347,-121.924252,29.3,1,1
H05247,36.896642,-121.922,29.3,1,1
H05247,36.896964,-121.93055,27.4,1,1
H05247,36.897017,-121.932008,29.3,1,1
H05247,36.897422,-121.92818,29.3,1,1
H05247,36.897439,-121.920791,25.6,1,1
H05247,36.897945,-121.924013,25.6,1,1
H05247,36.898231,-121.933338,31.1,1,1
H05247,36.898478,-121.930527,27.4,1,1


Get this into R and convert to SpatialPointsDataFrame

> survey = read.csv("./Data/H05247.xyz")
> coordinates(survey) = ~long + lat
> proj4string(survey) = CRS("+init=epsg:4326")
> survey = spTransform(survey, osm())
> plot(monterey)
> points(survey, pch = ".")
plot of chunk unnamed-chunk-10

Kriging

Interpolation of sampled values onto a grid.

> require(automap)
> depth = autoKrige(depth ~ 1, survey)
[using ordinary kriging]
> plot(depth)
plot of chunk unnamed-chunk-12

Mapping the output

Convert to a multi-layer raster object...

> depthR = stack(depth$krige_output)
> depthR
class       : RasterStack 
dimensions  : 77, 84, 6468, 3  (nrow, ncol, ncell, nlayers)
resolution  : 269.5, 269.5  (x, y)
extent      : -13583183, -13560542, 4408833, 4429588  (xmin, xmax, ymin, ymax)
coord. ref. : +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 
names       : var1.pred, var1.var, var1.stdev 
min values  :     14.89,   161.47,      12.71 
max values  :     744.4,    610.2,       24.7 
> plot(depthR[["var1.pred"]])
> points(survey, pch = 19, cex = 0.1)
plot of chunk unnamed-chunk-13

Standard Errors

> plot(depthR[["var1.stdev"]])
plot of chunk unnamed-chunk-14

Raster operations

> pR = raster(depthR[["var1.pred"]])
> pR[] = 1 - pnorm(200, values(depthR[["var1.pred"]]), 
     values(depthR[["var1.stdev"]]))
> plot(pR)
plot of chunk unnamed-chunk-15

Plot area of uncertainty:

> plot((pR < 0.999) & (pR > 0.001))
> contour(depthR[["var1.pred"]], 
     add = TRUE)
plot of chunk unnamed-chunk-16

Another survey

There is another survey data file (H05405) - repeat the analysis!

plot of chunk unnamed-chunk-18