1 References

3 Coordinate Reference Systems (CRS)

  • unprojected or geographic: Latitude and Longitude for referencing location on the ellipsoid Earth
    (decimal degrees (DD) or degrees, minutes, and seconds (DMS))

  • projected: Easting and Northing for referencing location on 2-dimensional representation of Earth
    Common projection: Universal Transverse Mercator (UTM)
    Location is given by the zone number (60 zones), hemisphere (north or south), and Easting and Northing coordinates in the zone in meters

4 Spatial data

Spatial data can be represented as simple features (sf) or as spatial class (sp).

Packages sf and sp can be used to work with spatial data.

4.1 Simple features (sf)

sf objects are extensions of data frames. Each row is one feature or spatial object that could have associated data.

4.2 Spatial class (sp)

  • SpatialPoints and SpatialPointsDataFrame
  • SpatialLines and SpatialLinesDataFrame
  • SpatialPolygons and SpatialPolygonsDataFrame
  • SpatialGrid and SpatialGridDataFrame
  • SpatialPixel and SpatialPixelDataFrame

5 Shapefile

Geographic data can be represented using a data storage format called shapefile. A shapefile consists of a collection of related files:

  • .shp: contains the geometry data,
  • .shx: is a positional index of the geometry data that allows to seek forwards and backwards the .shp file,
  • .dbf: stores the attributes for each shape.

Other files that can form a shapefile are the following:

  • .prj: plain text file describing the projection,
  • .sbn and .sbx: spatial index of the geometry data,
  • .shp.xml: geospatial metadata in XML format.

We can read shapefiles with rgdal::readOGR() or sf::st_read()

# name of the shapefile of North Carolina of the sf package
nameshp <- system.file("shape/nc.shp", package = "sf")


Exercise: Identify the folder of nc.shp in your computer.

5.1 sf (simple features)

# read shapefile with st_read()
library(sf)
## Warning: package 'sf' was built under R version 4.0.5
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
map <- st_read(nameshp, quiet = TRUE)
class(map)
## [1] "sf"         "data.frame"
head(map)
plot(map)
## Warning: plotting the first 10 out of 14 attributes; use max.plot = 14 to plot
## all

5.2 sp

# read shapefile with readOGR()
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/MORAGAPE/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/MORAGAPE/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/MORAGAPE/Documents/R/win-library/4.0/rgdal/proj
map <- rgdal::readOGR(nameshp, verbose = FALSE)
class(map)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
head(map@data)
plot(map)

spplot(map, zcol = "BIR74")

6 Types of spatial data

https://www.paulamoraga.com/book-geospatial/sec-spatialdataandCRS.html#types-of-spatial-data

Exercises:

  • Change theme
  • Add labels
  • Save plot

7 Making maps with R

https://www.paulamoraga.com/book-geospatial/sec-spatialdataandCRS.html#making-maps-with-r

Exercises:

Create map with leaflet or mapview depicting locations

library(leaflet)
## Warning: package 'leaflet' was built under R version 4.0.5
leaflet() %>% addTiles() %>% addMarkers(lng = 174.768, lat = -36.852, popup = "Location")

11 Exercises

11.1 Download map Saudi Arabia

library(rgeoboundaries)
library(ggplot2)
library(viridis)
## Loading required package: viridisLite
library(mapview)
library(rnaturalearth)

# map <- rgeoboundaries::geoboundaries(country = "SAU")
map <- ne_countries(type = "countries", country = "Saudi Arabia", scale = "medium", returnclass = "sf")
plot(map)

# map1 <- geoboundaries("Saudi Arabia", "adm1")
map1 <- rnaturalearth::ne_states("Saudi Arabia", returnclass = "sf")
names(map1)
##  [1] "featurecla" "scalerank"  "adm1_code"  "diss_me"    "iso_3166_2"
##  [6] "wikipedia"  "iso_a2"     "adm0_sr"    "name"       "name_alt"  
## [11] "name_local" "type"       "type_en"    "code_local" "code_hasc" 
## [16] "note"       "hasc_maybe" "region"     "region_cod" "provnum_ne"
## [21] "gadm_level" "check_me"   "datarank"   "abbrev"     "postal"    
## [26] "area_sqkm"  "sameascity" "labelrank"  "name_len"   "mapcolor9" 
## [31] "mapcolor13" "fips"       "fips_alt"   "woe_id"     "woe_label" 
## [36] "woe_name"   "latitude"   "longitude"  "sov_a3"     "adm0_a3"   
## [41] "adm0_label" "admin"      "geonunit"   "gu_a3"      "gn_id"     
## [46] "gn_name"    "gns_id"     "gns_name"   "gn_level"   "gn_region" 
## [51] "gn_a1_code" "region_sub" "sub_code"   "gns_level"  "gns_lang"  
## [56] "gns_adm1"   "gns_region" "min_label"  "max_label"  "min_zoom"  
## [61] "wikidataid" "name_ar"    "name_bn"    "name_de"    "name_en"   
## [66] "name_es"    "name_fr"    "name_el"    "name_hi"    "name_hu"   
## [71] "name_id"    "name_it"    "name_ja"    "name_ko"    "name_nl"   
## [76] "name_pl"    "name_pt"    "name_ru"    "name_sv"    "name_tr"   
## [81] "name_vi"    "name_zh"    "ne_id"      "geometry"
# print(map1, n = Inf)

map1$ISO <- map1$iso_3166_2

map1$vble <- rnorm(nrow(map1), 0, 10)

ggplot(map1) + geom_sf(aes(fill = vble)) + scale_fill_viridis() + theme_bw()

ggplot(map1) + geom_sf(aes(fill = ISO)) + scale_fill_viridis(discrete = TRUE) + theme_bw()

mapview(map1, zcol = "ISO", layer.name = "ISO")
mapview(map1, zcol = c("vble", "ISO"), layer.name = "ISO")

12 Retrieve spatial data

library(raster)
## Warning: package 'raster' was built under R version 4.0.5
## 
## Attaching package: 'raster'
## The following object is masked from 'package:MASS':
## 
##     select
temp <- getData(name = "worldclim", var = "tmean", res = 10)
temp <- temp*.1
plot(temp)

temp <- mean(temp)
plot(temp)

# Crop to region boundary
r2 <- crop(temp, extent(map1))
plot(r2)

r3 <- mask(r2, map1)
plot(r3)

12.1 Bounding box

#plot(sf::as_Spatial(map))
#axis(1)
#axis(2)

bb <- st_bbox(map)
bb
##     xmin     ymin     xmax     ymax 
## 34.61621 16.37178 55.64102 32.12451

12.2 Grid points within bounding box

# Grid of points encompassing the bounding box
x <- seq(bb$xmin, bb$xmax, length.out = 50)
y <- seq(bb$ymin, bb$ymax, length.out = 50)
coop <- as.matrix(expand.grid(x, y))
plot(coop)

12.3 Plot points with mapview

library(mapview)
coopsf <- st_as_sf(as.data.frame(coop), coords = c(1, 2), crs = 4326)
mapview(list(map, coopsf), layer.name = c("map", "points"), layers.control.pos = "topright")
## Warning in cbind(`Feature ID` = fid, mat): number of rows of result is not a
## multiple of vector length (arg 1)

12.4 Identifying points within map

coopsf$inside <- ifelse(sf::st_intersects(coopsf, map), 1, 0)
coopsf <- coopsf[which(coopsf$inside == 1), ]

# Plot with mapview
mapview(list(map, coopsf), layer.name = c("map", "points"), layers.control.pos = "topright")