Background

Analysis of Spatial Data in R

R has a rich package ecosystem for read/writing, manipulating, and analyzing spatial data.


Some core packages:

  • sp - core classes for handling spatial data.

  • rgdal - R interface to gdal (Geospatial Data Abstraction Library) for reading and writing spatial data.

  • maptools - Additional tools for reading and writing spatial data.

  • rgeos - R interface to geos (Geometry Engine Open Source) library for querying and manipulating spatial data. Reading and writing WKT.

  • raster - classes and tools for handling spatial raster data.


See more - Spatial task view

GeoJSON

Point

s = '{
        "type": "Point", 
        "coordinates": [30, 10]
     }'
sp = readOGR(s, "OGRGeoJSON", verbose = FALSE)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
points(coordinates(sp), pch=16)

plot of chunk unnamed-chunk-2

MultiPoint

s = '{
        "type": "MultiPoint", 
        "coordinates": [ [10, 40], [40, 30], [20, 20], [30, 10] ]
     }'
sp = readOGR(s, "OGRGeoJSON", verbose = FALSE)
## Warning: eType not chosen
## Error: Incompatible geometry: 4

LineString

s = '{
        "type": "LineString", 
        "coordinates": [ [30, 10], [10, 30], [40, 40] ]
     }'
sp = readOGR(s, "OGRGeoJSON", verbose = FALSE)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
r=rapply(coordinates(sp), points, pch=16)

plot of chunk unnamed-chunk-4

MultiLineString

s = '{
        "type": "MultiLineString", 
        "coordinates": [ [[10, 10], [20, 20], [10, 40]], 
                         [[40, 40], [30, 30], [40, 20], [30, 10]]
                       ]
     }'
sp = readOGR(s, "OGRGeoJSON", verbose = FALSE)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
r=rapply(coordinates(sp), points, pch=16)

plot of chunk unnamed-chunk-5

Polygon

s = '{
        "type": "Polygon", 
        "coordinates": [ [[30, 10], [40, 40], [20, 40], [10, 20], [30, 10]] ]
     }'
sp = readOGR(s, "OGRGeoJSON", verbose = FALSE)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-7

Polygon with Hole

s = '{
        "type": "Polygon", 
        "coordinates": [ [[35, 10], [45, 45], [15, 40], [10, 20], [35, 10]], 
                         [[20, 30], [35, 35], [30, 20], [20, 30]]
                       ]
     }'
sp = readOGR(s, "OGRGeoJSON", verbose = FALSE)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-8

MultiPolygon

s = '{
        "type": "MultiPolygon", 
        "coordinates": [ [ [[30, 20], [45, 40], [10, 40], [30, 20]] ], 
                         [ [[15, 5], [40, 10], [10, 20], [5, 10], [15, 5]] ]
                       ]
     }'
sp = readOGR(s, "OGRGeoJSON", verbose = FALSE)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-9

MultiPolygon with Hole(s)

s = '{
        "type": "MultiPolygon", 
        "coordinates": [  [ [[40, 40], [20, 45], [45, 30], [40, 40]] ], 
                          [ [[20, 35], [10, 30], [10, 10], [30, 5], [45, 20], [20, 35]], 
                            [[30, 20], [20, 15], [20, 25], [30, 20]]
                          ]
                       ]
     }'
sp = readOGR(s, "OGRGeoJSON", verbose = FALSE)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-10

GeometryCollection

s = '{
        "type": "GeometryCollection",
        "geometries": [
          { 
            "type": "Point",
            "coordinates": [30, 10]
          },
          { 
            "type": "LineString",
            "coordinates": [ [30, 10], [10, 30], [40, 40] ]
          }
        ]
     }'

sp = readOGR(s, "OGRGeoJSON", verbose = FALSE)
## Warning: eType not chosen
## Error: Incompatible geometry: 7

Well Known Text

Point

s = 'POINT (30 10)'
sp = readWKT(s)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
points(coordinates(sp), pch=16)

plot of chunk unnamed-chunk-13

MultiPoint

s = 'MULTIPOINT ((10 40), (40 30), (20 20), (30 10))'
sp = readWKT(s)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
points(coordinates(sp), pch=16)

plot of chunk unnamed-chunk-14

LineString

s = 'LINESTRING (30 10, 10 30, 40 40)'
sp = readWKT(s)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
r=rapply(coordinates(sp), points, pch=16)

plot of chunk unnamed-chunk-15

MultiLineString

s = 'MULTILINESTRING ((10 10, 20 20, 10 40), 
                      (40 40, 30 30, 40 20, 30 10))'
sp = readWKT(s)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
r=rapply(coordinates(sp), points, pch=16)

plot of chunk unnamed-chunk-16

Polygon

s = 'POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))'
sp = readWKT(s)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-17

Polygon with Hole

s = 'POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10),
              (20 30, 35 35, 30 20, 20 30))'
sp = readWKT(s)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-18

MultiPolygon

s = 'MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)),
                   ((15 5, 40 10, 10 20, 5 10, 15 5)))'
sp = readWKT(s)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-19

MultiPolygon with Hole(s)

s = 'MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)),
                   ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35),
                    (30 20, 20 15, 20 25, 30 20)))'
sp = readWKT(s)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-20

GeometryCollection

s = 'GEOMETRYCOLLECTION (POINT (4 8),
                         LINESTRING (4 6,7 10),
                         POLYGON ((6 6, 8 6, 8 8, 6 6)))'

sp = readWKT(s)

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)

plot of chunk unnamed-chunk-21

sp

Point

sp = SpatialPoints(data.frame(x=30,y=10))

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
points(coordinates(sp), pch=16)

plot of chunk unnamed-chunk-22

MultiPoint

sp = SpatialPoints(data.frame(x=c(10,40,20,30),y=c(40,30,20,10)))

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
points(coordinates(sp), pch=16)

plot of chunk unnamed-chunk-23

LineString

l = Line(data.frame(x=c(30,10,40), y=c(10,30,40)))

ls = Lines(list(l), ID=1)
sp = SpatialLines(list(ls))

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
r=rapply(coordinates(sp), points, pch=16)

plot of chunk unnamed-chunk-24

MultiLineString

l1 = Line(data.frame(x=c(10,20,10), y=c(10,20,40)))
l2 = Line(data.frame(x=c(40,30,40,30), y=c(40,30,20,10)))

ls = Lines(list(l1,l2), ID=1)
sp = SpatialLines(list(ls))

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE)
r=rapply(coordinates(sp), points, pch=16)

plot of chunk unnamed-chunk-25

Polygon

p = Polygon(data.frame(x=c(30,40,20,10,30), y=c(10,40,40,20,10)))

ps = Polygons(list(p), ID=1)
sp = SpatialPolygons(list(ps))

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-26

Polygon with Hole

p = Polygon(data.frame(x=c(35,45,15,10,35), y=c(10,45,40,20,10)))
h = Polygon(data.frame(x=c(20,35,30,20), y=c(30,35,20,30)), hole=TRUE)

ps = Polygons(list(p,h), ID=1)
sp = SpatialPolygons(list(ps))

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-27

MultiPolygon

p1 = Polygon(data.frame(x=c(30,45,10,30), y=c(20,40,40,20)))
p2 = Polygon(data.frame(x=c(15,40,10,5,15), y=c(5,10,20,10,5)))

ps = Polygons(list(p1,p2), ID=1)
sp = SpatialPolygons(list(ps))

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-28

MultiPolygon with Hole(s)

p1 = Polygon(data.frame(x=c(40,20,45,40), y=c(40,45,30,40)))
p2 = Polygon(data.frame(x=c(20,10,10,30,45,20), y=c(35,30,10,5,20,35)))
h1 = Polygon(data.frame(x=c(30,20,20,30), y=c(20,15,25,20)), hole=TRUE)

ps = Polygons(list(p1,p2,h1), ID=1)
sp = SpatialPolygons(list(ps))

par(mar=c(4,4,1,1))
plot(sp, axes=TRUE, col="lightgrey")
points(poly_coords(sp), pch=16)

plot of chunk unnamed-chunk-29

Exercise

Largest Island in a lake on an island in a lake on an island