Geographic data corresponds to data that have been linked to a geographic coordinate system
Working with geographically located data is a relatively new phenomenon
Nowadays such lack of geographic data is hard to imagine
## install.packages("leaflet")
library(leaflet)
popup = c("Storrs, Connecticut")
leaflet() %>%
setView(lng = -72.251194, lat = 41.807417, zoom = 1)%>%
addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") %>%
addMarkers(lng = c(-72.251194),
lat = c(41.807417),
popup = popup)popup = c("SSH Building")
leaflet() %>%
setView(lng = -72.251194, lat = 41.807417, zoom = 12)%>%
addTiles() %>%
addMarkers(lng = c(-72.251194),
lat = c(41.807417),
popup = popup)There are many ways to handle geographic data in R, with dozens of packages in the area. The ones that we will focus on in this course are:
spsf (very functional with the tidyverse)rasterHonorable mentions:
maptoolsrgdalsp PackageThe sp package provides the original infrastructure for handling spatial data in R. It defines formal S4 classes for representing geographic features—points, lines, polygons, and grids—together with their coordinate reference systems (CRS) and attribute data. These structures form the foundation for spatial analysis, visualization, and conversion between coordinate systems.
Although newer packages like sf now offer simpler data-frame-based representations, sp remains widely used and is essential for understanding many legacy R spatial workflows.
sp.sp ObjectLet’s build a sp object:
#install.packages("sp")
library(sp)
#create the highways
ln1 <- Line(matrix(runif(6), ncol=2))
ln1An object of class "Line"
Slot "coords":
[,1] [,2]
[1,] 0.9089826 0.1298057
[2,] 0.2503783 0.4034792
[3,] 0.9638094 0.4657451ln2 <- Line(matrix(runif(6), ncol=2))
ln2An object of class "Line"
Slot "coords":
[,1] [,2]
[1,] 0.4099365 0.6904562
[2,] 0.7023185 0.9883888
[3,] 0.8222399 0.9577005str(ln1)Formal class 'Line' [package "sp"] with 1 slot
..@ coords: num [1:3, 1:2] 0.909 0.25 0.964 0.13 0.403 ...Note that the @coords slot holds the coordinates
## combine the highways to a Lines object
lns1 <- Lines(list(ln1), ID = c("hwy1"))
lns2 <- Lines(list(ln2), ID = c("hwy2"))
str(lns1)Formal class 'Lines' [package "sp"] with 2 slots
..@ Lines:List of 1
.. ..$ :Formal class 'Line' [package "sp"] with 1 slot
.. .. .. ..@ coords: num [1:3, 1:2] 0.909 0.25 0.964 0.13 0.403 ...
..@ ID : chr "hwy1"The Line objects are now in a list and we have an additional ID slot, which uniquely identifies each Line object .
# create a geospatial object with SpatialLines
sp_lns <- SpatialLines(list(lns1, lns2))
str(sp_lns)Formal class 'SpatialLines' [package "sp"] with 3 slots
..@ lines :List of 2
.. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
.. .. .. ..@ Lines:List of 1
.. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
.. .. .. .. .. .. ..@ coords: num [1:3, 1:2] 0.909 0.25 0.964 0.13 0.403 ...
.. .. .. ..@ ID : chr "hwy1"
.. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
.. .. .. ..@ Lines:List of 1
.. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
.. .. .. .. .. .. ..@ coords: num [1:3, 1:2] 0.41 0.702 0.822 0.69 0.988 ...
.. .. .. ..@ ID : chr "hwy2"
..@ bbox : num [1:2, 1:2] 0.25 0.13 0.964 0.988
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "x" "y"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr NAA data frame can be created with a column of SpatialLines.
## create a SpatialLinesDataframe
dfr <- data.frame(id = c("hwy1", "hwy2"), # note how we use the same IDs from above!
cars_per_hour = c(78, 22))
sp_lns_dfr <- SpatialLinesDataFrame(sp_lns, dfr, match.ID = "id")
str(sp_lns_dfr)Formal class 'SpatialLinesDataFrame' [package "sp"] with 4 slots
..@ data :'data.frame': 2 obs. of 2 variables:
.. ..$ id : chr [1:2] "hwy1" "hwy2"
.. ..$ cars_per_hour: num [1:2] 78 22
..@ lines :List of 2
.. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
.. .. .. ..@ Lines:List of 1
.. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
.. .. .. .. .. .. ..@ coords: num [1:3, 1:2] 0.909 0.25 0.964 0.13 0.403 ...
.. .. .. ..@ ID : chr "hwy1"
.. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
.. .. .. ..@ Lines:List of 1
.. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
.. .. .. .. .. .. ..@ coords: num [1:3, 1:2] 0.41 0.702 0.822 0.69 0.988 ...
.. .. .. ..@ ID : chr "hwy2"
..@ bbox : num [1:2, 1:2] 0.25 0.13 0.964 0.988
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "x" "y"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr NAspplot(sp_lns_dfr, "cars_per_hour")
spSpatial object — Any data object with geographic coordinates
(points, lines, polygons, or grids).
CRS (Coordinate Reference System) — Describes how coordinates map to Earth. Example:
CRS("+proj=longlat +datum=WGS84").
proj4string — Text form of CRS stored in an sp object.
@data — Attribute table (like columns in a data frame).
Access with obj@data.
@coords — Coordinate matrix of locations.
bbox — Bounding box (spatial extent) of the object.
spTransform() — Converts coordinates to another CRS.
over() — Spatial join (find which polygon a point belongs to).
Tip: The modern package
sfreplacesspand uses a simpler data frame structure with ageometrycolumn.
sf Packagesf implements a formal standard called “Simple Features” that specifies a storage and access model of spatial geometries (point, line, polygon).sp approach.Let’s recreate the sp object from above in sf
## install.packages("sf")
library(sf)
## create the highways
lnstr_sfg1 <- st_linestring(matrix(runif(6), ncol=2))
lnstr_sfg1
lnstr_sfg2 <- st_linestring(matrix(runif(6), ncol=2))
lnstr_sfg2
class(lnstr_sfg1)[1] "XY" "LINESTRING" "sfg" We would then combine this into a simple feature collection.
#combine the highways into a simple feature
lnstr_sfc <- st_sfc(lnstr_sfg1, lnstr_sfg2) # just one feature here
lnstr_sfcGeometry set for 2 features
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 0.05486331 ymin: 0.08513958 xmax: 0.9232197 ymax: 0.7123036
CRS: NA#create our spatial data frame
lnstr_sf <- st_sf(dfr , lnstr_sfc)
lnstr_sfplot(lnstr_sf["cars_per_hour"])
ggplot(lnstr_sf) +
geom_sf(aes(color = cars_per_hour))
Useful methods in sf:
methods(class = "sf") [1] [ [[<-
[3] [<- $<-
[5] aggregate anti_join
[7] arrange as.data.frame
[9] cbind coerce
[11] dbDataType dbWriteTable
[13] distinct dplyr_reconstruct
[15] drop_na duplicated
[17] filter full_join
[19] gather group_by
[21] group_split identify
[23] initialize inner_join
[25] left_join merge
[27] mutate nest
[29] pivot_longer pivot_wider
[31] plot points
[33] print rbind
[35] rename_with rename
[37] right_join rowwise
[39] sample_frac sample_n
[41] select semi_join
[43] separate_rows separate
[45] show slice
[47] slotsFromS3 spread
[49] st_agr st_agr<-
[51] st_area st_as_s2
[53] st_as_sf st_as_sfc
[55] st_bbox st_boundary
[57] st_break_antimeridian st_buffer
[59] st_cast st_centroid
[61] st_collection_extract st_concave_hull
[63] st_convex_hull st_coordinates
[65] st_crop st_crs
[67] st_crs<- st_difference
[69] st_drop_geometry st_exterior_ring
[71] st_filter st_geometry
[73] st_geometry<- st_inscribed_circle
[75] st_interpolate_aw st_intersection
[77] st_intersects st_is_full
[79] st_is_valid st_is
[81] st_join st_line_merge
[83] st_m_range st_make_valid
[85] st_minimum_bounding_circle st_minimum_rotated_rectangle
[87] st_nearest_points st_node
[89] st_normalize st_point_on_surface
[91] st_polygonize st_precision
[93] st_reverse st_sample
[95] st_segmentize st_set_precision
[97] st_shift_longitude st_simplify
[99] st_snap st_sym_difference
[101] st_transform st_triangulate_constrained
[103] st_triangulate st_union
[105] st_voronoi st_wrap_dateline
[107] st_write st_z_range
[109] st_zm summarise
[111] text transform
[113] transmute ungroup
[115] unite unnest
see '?methods' for accessing help and source codeRaster data represent spatial information as a grid of cells (or pixels), where each cell stores a value for a specific variable such as elevation, temperature, or land cover. The cells are arranged in rows and columns, covering a continuous area of space.
Raster data are typically used for continuous phenomena that vary smoothly across space. The spatial resolution depends on the cell size: smaller cells capture finer detail but require more storage.
Common raster formats include GeoTIFF and .grd. In R, raster data are handled by packages such as raster and terra.
## install.packages("raster")
library(raster)
HARV <- raster("data/HARV_RGB_Ortho.tif")
HARVclass : RasterLayer
band : 1 (of 3 bands)
dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
resolution : 0.25, 0.25 (x, y)
extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
source : HARV_RGB_Ortho.tif
names : HARV_RGB_Ortho_1
values : 0, 255 (min, max)plot(HARV)
This raster is generated as part of the NEON Harvard Forest field site .
## useful commands
crs(HARV)Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
WKT2 2019 representation:
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 18N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-75,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16018]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]] methods(class = class(HARV)) [1] ! != [
[4] [[ [[<- [<-
[7] %in% == $
[10] $<- addLayer adjacent
[13] aggregate all.equal area
[16] Arith as.array as.character
[19] as.data.frame as.factor as.integer
[22] as.list as.logical as.matrix
[25] as.raster as.vector asFactor
[28] atan2 bandnr barplot
[31] bbox blockSize boundaries
[34] boxplot brick buffer
[37] calc cellFromRowCol cellFromRowColCombine
[40] cellFromXY cellStats clamp
[43] click clump coerce
[46] colFromCell colFromX colSums
[49] Compare contour coordinates
[52] corLocal couldBeLonLat cover
[55] crop crosstab crs<-
[58] cut cv density
[61] dim dim<- direction
[64] disaggregate distance extend
[67] extent extract flip
[70] focal freq getValues
[73] getValuesBlock getValuesFocal gridDistance
[76] hasValues head hist
[79] image init initialize
[82] inMemory interpolate intersect
[85] is.factor is.finite is.infinite
[88] is.na is.nan isLonLat
[91] KML labels layerize
[94] length levels levels<-
[97] lines localFun log
[100] Logic mask match
[103] Math Math2 maxValue
[106] mean merge metadata
[109] minValue modal mosaic
[112] names names<- ncell
[115] ncol ncol<- nlayers
[118] nrow nrow<- origin
[121] origin<- overlay persp
[124] plot predict print
[127] proj4string proj4string<- quantile
[130] raster rasterize ratify
[133] readAll readStart readStop
[136] reclassify rectify res
[139] res<- resample RGB
[142] rotate rowColFromCell rowFromCell
[145] rowFromY rowSums sampleRandom
[148] sampleRegular sampleStratified scale
[151] select setMinMax setValues
[154] shift show spplot
[157] stack stackSelect stretch
[160] subs subset Summary
[163] summary t tail
[166] terrain text trim
[169] unique update values
[172] values<- Which which.max
[175] which.min wkt writeRaster
[178] writeStart writeStop writeValues
[181] xFromCell xFromCol xmax
[184] xmax<- xmin xmin<-
[187] xres xyFromCell yFromCell
[190] yFromRow ymax ymax<-
[193] ymin ymin<- yres
[196] zonal zoom
see '?methods' for accessing help and source codeWe are going to create a few different types of maps, which will require different R packages.
Static Maps
Most common type of visual output
plot() most commonly used command from base are
Can utilize ggplot() and tmap()
Animated Maps
Are becoming more widely used than in the past
tmap()Interactive Maps
While static and animated maps can enliven geographic datasets, interactive maps can take them to a new level
leaflet()sf)sf package extends the base plot command, so it can be used on sf objects.sflibrary(sf)
philly_crimes_sf <- st_read("data/PhillyCrimerate.shp")Reading layer `PhillyCrimerate' from data source
`/Users/junyan/work/teaching/1010-f25/1010f25/data/PhillyCrimerate.shp'
using driver `ESRI Shapefile'
Simple feature collection with 384 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 1739497 ymin: 457343.7 xmax: 1764030 ymax: 490544.9
Projected CRS: Albersglimpse(philly_crimes_sf)Rows: 384
Columns: 5
$ GEOID10 <chr> "42101000100", "42101000200", "42101000300", "42101000401",…
$ n_homic <int> 10, 4, 7, 9, 6, 9, 1, 3, 1, 1, 1, 2, 1, 1, 4, 1, 1, 1, 3, 2…
$ tract_area <dbl> 704688.5, 381197.4, 550949.2, 231327.5, 303868.9, 427352.7,…
$ homic_rate <dbl> 14.190667, 10.493252, 12.705346, 38.905882, 19.745355, 21.0…
$ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((1752148 467..., MULTIPOLYGON (…plot(philly_crimes_sf)
Notice that the plot function will try and provide a map of every variable in the data set.
To plot a single attribute, we need to provide an object of class sf.
plot(philly_crimes_sf$homic_rate) # this is a numeric vector!
plot(philly_crimes_sf["homic_rate"])
Based on the plot above, it would seem that we have a lot of “empty” geographic units or unevenly distributed
ggplot(philly_crimes_sf) +
geom_histogram(mapping = aes(x = homic_rate))
Let’s break the plot up into quantiles to better distinguish the census tracts with low values.
plot(philly_crimes_sf["homic_rate"],
main = "Philadelphia homicide density per square km",
breaks = "quantile")
ggplot()Package ggplot is a widely used and powerful plotting library for R as we have utilized over the course of the semester. It is not specifically geared towards mapping, but one can generate great maps.
ggplot(philly_crimes_sf) +
geom_sf(aes(fill=homic_rate))
Homicide rate is a continuous variable and is plotted by ggplot as such. If we wanted to plot our map as a ‘true’ choropleth map we need to convert our continuous variable into a categorical one, according to whichever brackets we want to use.
breaks_qt <- classInt::classIntervals(
c(min(philly_crimes_sf$homic_rate) - .00001,
philly_crimes_sf$homic_rate),
n = 7, style = "quantile")
philly_crimes_sf <- mutate(philly_crimes_sf,
homic_rate_cat = cut(homic_rate, breaks_qt$brks))
ggplot(philly_crimes_sf) +
geom_sf(aes(fill=homic_rate_cat)) +
scale_fill_brewer(palette = "OrRd") 
tmap()Package tmap is based on the idea of a ‘grammar of graphics’ similar to ggplot This involves a separation between the input data and the aesthetics (how data are visualized): each input dataset can be ‘mapped’ in a range of different ways including location on the map (defined by data’s geometry), color, and other visual variables. The package is specifically designed to make creation of thematic maps more convenient with ggplot syntax.
tm_shape() defines the data, raster, and vector objects.
tm_polygons() sets the attribute variable to map, the break style, and a title.
tm_fill() & tm_borders() control the layers.
library(tmap)
m <- tm_shape(philly_crimes_sf) +
tm_polygons(
"homic_rate",
fill.scale = tm_scale_intervals(style = "quantile", n = 4),
fill.legend = tm_legend("Philadelphia \nhomicide density \nper sqKm")
)
m
Package tmap has a nice feature that allows us to give basic interactivity to the map. We can switch from plot mode into view mode and call the plot.
tmap_mode("view")
mLayers can be controlled by tm_fill() & tm_borders().
#
#install.packages("spData")
library(spData)
library(tmap)
# Add fill layer to nz shape
tm_shape(nz) +
tm_fill() # Add border layer to nz shape
tm_shape(nz) +
tm_borders() # Add fill and border layers to nz shape
tm_shape(nz) +
tm_fill() +
tm_borders() A useful feature of tmap is its ability to store objects representing maps
map_nz = tm_shape(nz) + tm_polygons()
class(map_nz)[1] "tmap"map_nz1 <- map_nz +
tm_polygons("Population", fill.scale = tm_scale_intervals(),
palette = "viridis",
title = "Population (by region)") +
tm_layout(frame = FALSE)
map_nz1There are no limits to the number of layers or shapes that can be added to tmap objects. Additionally, tmap allows multiple map objects to be arranged in a single plot with tmap_arrange()
map_nz2 <- tm_shape(nz) +
tm_polygons("Median_income", style = "pretty", palette = "brewer.blues",
title = "Median income")
# combined
tmap_arrange(map_nz1, map_nz2) #, map_nz3)gifskiConsider the United Nations urban populaton data including projections. Package tmap makes it easy to create maps with year as facet.
urb_1970_2030 <- urban_agglomerations |>
filter(year == 1970| year == 1990 | year == 2010 | year == 2030)
tm_shape(world) +
tm_polygons() +
tm_shape(urb_1970_2030) +
tm_symbols(col = "black", border.col = "white", size = "population_millions") +
tm_facets(by = "year", nrow = 2, free.coords = FALSE)
urb_anim = tm_shape(world) + tm_polygons() +
tm_shape(urban_agglomerations) + tm_dots(size = "population_millions") +
tm_facets(along = "year", free.coords = FALSE)
Package gifski takes the urb_anim object and create a GIF image, in which each year is a frame.
library(gifski)
tmap_animation(urb_anim, filename = "urb_anim.gif", delay = 25)leaflet()Package leaflet provides bindings to the Leaflet JavaScript library, “the leading open-source JavaScript library for mobile-friendly interactive maps.” We have already seen a simple use of leaflet in the tmap example.
Let’s build up the map step by step
library(leaflet)
# reproject
philly_WGS84 <- st_transform(philly_crimes_sf, 4326)
leaflet(philly_WGS84) %>%
addPolygons()To map the homicide density we use addPolygons()
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)
p_popup <- paste0("<strong>Homicide Density: </strong>", philly_WGS84$homic_rate)
m <- leaflet(philly_WGS84) %>%
addPolygons(
stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun(homic_rate), # set fill color with function from above and value
fillOpacity = 0.8, smoothFactor = 0.5, # make it nicer
popup = p_popup) # add popup
mWe can add a basemap, which defaults to OSM, with addTiles()
m <- addTiles(m)
mLastly, we add a legend with addLegend()
m <- addLegend(m,
"bottomright", # location
pal=pal_fun, # palette function
values=~homic_rate, # value to be passed to palette function
title = 'Philadelphia homicide density per sqkm') # legend title
m