Spatial Data

This document introduces various R packages for representing spatial data in R. It also explores packages for reading/writing spatial data from external formats.

Vecor Data Using sp Package

Construct a sp object.

Read the meuse data.frame and convert it into a spatial object. Explore the result.

library(sp)

Loading required package: methods

data(meuse)
print(class(meuse))

[1] "data.frame"

print(colnames(meuse))

 [1] "x"       "y"       "cadmium" "copper"  "lead"    "zinc"    "elev"   
 [8] "dist"    "om"      "ffreq"   "soil"    "lime"    "landuse" "dist.m" 

{{ coordinates(meuse) <- ~x+y }}
{{ proj4string(meuse) <- CRS("+init=epsg:28992") }}
print(class(meuse))

[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

summary(meuse)

Object of class SpatialPointsDataFrame
Coordinates:
     min    max
x 178605 181390
y 329714 333611
Is projected: TRUE 
proj4string :
[+init=epsg:28992 +proj=sterea +lat_0=52.15616055555555
+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000
+ellps=bessel
+towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725
+units=m +no_defs]
Number of points: 155
Data attributes:
    cadmium           copper            lead            zinc       
 Min.   : 0.200   Min.   : 14.00   Min.   : 37.0   Min.   : 113.0  
 1st Qu.: 0.800   1st Qu.: 23.00   1st Qu.: 72.5   1st Qu.: 198.0  
 Median : 2.100   Median : 31.00   Median :123.0   Median : 326.0  
 Mean   : 3.246   Mean   : 40.32   Mean   :153.4   Mean   : 469.7  
 3rd Qu.: 3.850   3rd Qu.: 49.50   3rd Qu.:207.0   3rd Qu.: 674.5  
 Max.   :18.100   Max.   :128.00   Max.   :654.0   Max.   :1839.0  
                                                                   
      elev             dist               om         ffreq  soil   lime   
 Min.   : 5.180   Min.   :0.00000   Min.   : 1.000   1:84   1:97   0:111  
 1st Qu.: 7.546   1st Qu.:0.07569   1st Qu.: 5.300   2:48   2:46   1: 44  
 Median : 8.180   Median :0.21184   Median : 6.900   3:23   3:12          
 Mean   : 8.165   Mean   :0.24002   Mean   : 7.478                        
 3rd Qu.: 8.955   3rd Qu.:0.36407   3rd Qu.: 9.000                        
 Max.   :10.520   Max.   :0.88039   Max.   :17.000                        
                                    NA's   :2                             
    landuse       dist.m      
 W      :50   Min.   :  10.0  
 Ah     :39   1st Qu.:  80.0  
 Am     :22   Median : 270.0  
 Fw     :10   Mean   : 290.3  
 Ab     : 8   3rd Qu.: 450.0  
 (Other):25   Max.   :1000.0  
 NA's   : 1                   

Load sp objects from packages.

Load all ‘meuse’ datasets

demo("meuse", ask=FALSE)

    demo(meuse)
    ---- ~~~~~


✔  user2017.geodataviz master*
❯ require(sp)


✔  user2017.geodataviz master*
❯ crs = CRS("+init=epsg:28992")


✔  user2017.geodataviz master*
❯ data("meuse")


✔  user2017.geodataviz master*
❯ coordinates(meuse) <- ~x+y


✔  user2017.geodataviz master*
❯ proj4string(meuse) <- crs


✔  user2017.geodataviz master*
❯ data("meuse.grid")


✔  user2017.geodataviz master*
❯ coordinates(meuse.grid) <- ~x+y


✔  user2017.geodataviz master*
❯ gridded(meuse.grid) <- TRUE


✔  user2017.geodataviz master*
❯ proj4string(meuse.grid) <- crs


✔  user2017.geodataviz master*
❯ data("meuse.riv")


✔  user2017.geodataviz master*
❯ meuse.riv <- SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)),"meuse.riv")))


✔  user2017.geodataviz master*
❯ proj4string(meuse.riv) <- crs


✔  user2017.geodataviz master*
❯ data("meuse.area")


✔  user2017.geodataviz master*
❯ meuse.area = SpatialPolygons(list(Polygons(list(Polygon(meuse.area)), "area")))


✔  user2017.geodataviz master*
❯ proj4string(meuse.area) <- crs

Read sp objects from files

library(rgdal)
dsn <- system.file("vectors", package = "rgdal")[1]

ogrInfo(dsn=dsn, layer="scot_BNG")

Source: "/home/bhaskar/Library/R/3.x/cran/rgdal/vectors", layer: "scot_BNG"
Driver: ESRI Shapefile; number of rows: 56 
Feature type: wkbPolygon with 2 dimensions
Extent: (7094.552 529495) - (468285.5 1218342)
CRS: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs  
LDID: 0 
Number of fields: 13 
     name type length typeName
1   SP_ID    4      5   String
2    NAME    4     13   String
3    ID_x    2     19     Real
4   COUNT    2     19     Real
5     SMR    2     19     Real
6    LONG    2     19     Real
7     LAT    2     19     Real
8      PY    2     19     Real
9    EXP_    2     19     Real
10    AFF    2     19     Real
11 X_COOR    2     19     Real
12 Y_COOR    2     19     Real
13   ID_y    2     19     Real

print(OGRSpatialRef(dsn=dsn, layer="scot_BNG"))

[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs "

scot_BNG <- readOGR(dsn=dsn, layer="scot_BNG")

OGR data source with driver: ESRI Shapefile 
Source: "/home/bhaskar/Library/R/3.x/cran/rgdal/vectors", layer: "scot_BNG"
with 56 features
It has 13 fields

print(class(scot_BNG))

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

summary(scot_BNG)

Object of class SpatialPolygonsDataFrame
Coordinates:
         min       max
x   7094.552  468285.5
y 529495.039 1218342.5
Is projected: TRUE 
proj4string :
[+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
+y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy
+towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894]
Data attributes:
     SP_ID              NAME         ID_x           COUNT       
 0      : 1   Aberdeen    : 1   Min.   : 1.00   Min.   : 0.000  
 1      : 1   Angus       : 1   1st Qu.:14.75   1st Qu.: 4.750  
 10     : 1   Annandale   : 1   Median :28.50   Median : 8.000  
 11     : 1   Argyll-Bute : 1   Mean   :28.50   Mean   : 9.571  
 12     : 1   Banff-Buchan: 1   3rd Qu.:42.25   3rd Qu.:11.000  
 13     : 1   Bearsden    : 1   Max.   :56.00   Max.   :39.000  
 (Other):50   (Other)     :50                                   
      SMR             LONG            LAT              PY         
 Min.   :  0.0   Min.   :54.94   Min.   :1.430   Min.   :  27075  
 1st Qu.: 49.6   1st Qu.:55.78   1st Qu.:3.288   1st Qu.: 100559  
 Median :111.5   Median :56.04   Median :4.090   Median : 182333  
 Mean   :152.6   Mean   :56.40   Mean   :4.012   Mean   : 267498  
 3rd Qu.:223.0   3rd Qu.:57.02   3rd Qu.:4.730   3rd Qu.: 313845  
 Max.   :652.2   Max.   :60.24   Max.   :6.800   Max.   :2316353  
                                                                  
      EXP_             AFF             X_COOR           Y_COOR       
 Min.   : 1.100   Min.   : 0.000   Min.   :112892   Min.   : 561163  
 1st Qu.: 4.050   1st Qu.: 1.000   1st Qu.:256624   1st Qu.: 649520  
 Median : 6.300   Median : 7.000   Median :287577   Median : 681524  
 Mean   : 9.575   Mean   : 8.661   Mean   :288524   Mean   : 723127  
 3rd Qu.:10.125   3rd Qu.:11.500   3rd Qu.:333948   3rd Qu.: 794380  
 Max.   :88.700   Max.   :24.000   Max.   :442244   Max.   :1168904  
                                                                     
      ID_y      
 Min.   : 1.00  
 1st Qu.:14.75  
 Median :28.50  
 Mean   :28.50  
 3rd Qu.:42.25  
 Max.   :56.00  
                

Selecting/Subsetting

# 9th and 10th rows, and lead and zinc columns
meuse[9:10, c('lead','zinc')]

        coordinates lead zinc
9  (181060, 333231)  133  347
10 (181232, 333168)   80  183

# Filtering on a criteria
dim(meuse[meuse$lead>quantile(meuse$lead,.75),])

[1] 38 12

# Attribute data
range(meuse[['zinc']])

[1]  113 1839

Projection Transformation

print(proj4string(meuse))

[1] "+init=epsg:28992 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs"

meuse.latlon <- spTransform(meuse, CRS("+init=epsg:4326")) # What's '+init=epsg:4326' ?
print(proj4string(meuse))

[1] "+init=epsg:28992 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs"

Other Operations

• Vignettes 1 and 2 for sp.
• Packages gdistance, geosphere.

Vector Data Using sf Package

sf is the successor to sp and provides much faster and tidyverse compatible API.

Read a sf Object

library(sf)

Linking to GEOS 3.5.1, GDAL 2.1.3, proj.4 4.9.2, lwgeom 2.3.2 r15302

scot_BNG <- sf::st_read(system.file("vectors", package = "rgdal")[1], 'scot_BNG')

Reading layer `scot_BNG' from data source `/home/bhaskar/Library/R/3.x/cran/rgdal/vectors' using driver `ESRI Shapefile'
Simple feature collection with 56 features and 13 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 7094.552 ymin: 529495 xmax: 468285.5 ymax: 1218342
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs

print(class(scot_BNG))

[1] "sf"         "data.frame"

summary(scot_BNG)

     SP_ID              NAME         ID_x           COUNT       
 0      : 1   Aberdeen    : 1   Min.   : 1.00   Min.   : 0.000  
 1      : 1   Angus       : 1   1st Qu.:14.75   1st Qu.: 4.750  
 10     : 1   Annandale   : 1   Median :28.50   Median : 8.000  
 11     : 1   Argyll-Bute : 1   Mean   :28.50   Mean   : 9.571  
 12     : 1   Banff-Buchan: 1   3rd Qu.:42.25   3rd Qu.:11.000  
 13     : 1   Bearsden    : 1   Max.   :56.00   Max.   :39.000  
 (Other):50   (Other)     :50                                   
      SMR             LONG            LAT              PY         
 Min.   :  0.0   Min.   :54.94   Min.   :1.430   Min.   :  27075  
 1st Qu.: 49.6   1st Qu.:55.78   1st Qu.:3.288   1st Qu.: 100559  
 Median :111.5   Median :56.04   Median :4.090   Median : 182333  
 Mean   :152.6   Mean   :56.40   Mean   :4.012   Mean   : 267498  
 3rd Qu.:223.0   3rd Qu.:57.02   3rd Qu.:4.730   3rd Qu.: 313845  
 Max.   :652.2   Max.   :60.24   Max.   :6.800   Max.   :2316353  
                                                                  
      EXP_             AFF             X_COOR           Y_COOR       
 Min.   : 1.100   Min.   : 0.000   Min.   :112892   Min.   : 561163  
 1st Qu.: 4.050   1st Qu.: 1.000   1st Qu.:256624   1st Qu.: 649520  
 Median : 6.300   Median : 7.000   Median :287577   Median : 681524  
 Mean   : 9.575   Mean   : 8.661   Mean   :288524   Mean   : 723127  
 3rd Qu.:10.125   3rd Qu.:11.500   3rd Qu.:333948   3rd Qu.: 794380  
 Max.   :88.700   Max.   :24.000   Max.   :442244   Max.   :1168904  
                                                                     
      ID_y                geometry 
 Min.   : 1.00   MULTIPOLYGON :56  
 1st Qu.:14.75   epsg:NA      : 0  
 Median :28.50   +proj=tmer...: 0  
 Mean   :28.50                     
 3rd Qu.:42.25                     
 Max.   :56.00                     
                                   

sf operations

# Just the geometry
scot_BNG.geom <- st_geometry(scot_BNG)
print(class(scot_BNG.geom))

[1] "sfc_MULTIPOLYGON" "sfc"             

# Reprojection
st_crs(scot_BNG)

$epsg
[1] NA

$proj4string
[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs"

attr(,"class")
[1] "crs"

scot_BNG.latlon <- st_transform(scot_BNG, '+init=epsg:4326')
st_crs(scot_BNG.latlon)

$epsg
[1] 4326

$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"

attr(,"class")
[1] "crs"

# Centroids
summary(scot_BNG$geometry)

 MULTIPOLYGON       epsg:NA +proj=tmer... 
           56             0             0 

scot_BNG.centers <- st_centroid(scot_BNG)
summary(scot_BNG.centers$geometry)

        POINT       epsg:NA +proj=tmer... 
           56             0             0 

Other Operations

methods(class='sf')

 [1] aggregate         cbind             coerce           
 [4] initialize        merge             plot             
 [7] print             rbind             [                
[10] [[<-              $<-               show             
[13] slotsFromS3       st_agr<-          st_agr           
[16] st_as_sf          st_bbox           st_boundary      
[19] st_buffer         st_cast           st_centroid      
[22] st_convex_hull    st_coordinates    st_crs<-         
[25] st_crs            st_difference     st_geometry<-    
[28] st_geometry       st_intersection   st_is            
[31] st_line_merge     st_make_valid     st_polygonize    
[34] st_precision      st_segmentize     st_set_precision 
[37] st_simplify       st_split          st_sym_difference
[40] st_transform      st_triangulate    st_union         
[43] st_voronoi        st_zm            
see '?methods' for accessing help and source code

Raster Data using raster Package

Read a Raster file

library(raster)

Attaching package: 'raster'

The following object is masked from 'package:magrittr':

    extract

r <- raster(system.file("external/test.grd", package="raster"))
print(class(r))

[1] "RasterLayer"
attr(,"package")
[1] "raster"

summary(r)

             test
Min.     128.4340
1st Qu.  293.2325
Median   371.4120
3rd Qu.  499.8195
Max.    1805.7800
NA's    6097.0000

raster Operations

r_log <- calc(r, log)
par(mar=c(5,5,5,5), mfrow=c(1,2))
hist(r); hist(r_log)

For many more see the raster package vignette.