This chapter introduces a number of concepts associated with handling spatial data, and points forward to later sections where they are discussed in more detail.
The code sections in this chapter are not explained, but can be unfolded by clicking on the “code” button. They should be easy to follow when understanding R at the level of, say, R for Data Science (Wickham and Grolemund 2017). Detailed explanation of R code starts in Part II of this book.
The typical way to graph spatial data is by creating a map. Let us consider a simple map, shown in figure 1.1.
A number of graphical elements are present here, in this case:
- polygons are drawn with a black outline and filled with colors chosen according to a variable
BIR74, whose name is in the title
- a legend key explains the meaning of the colors, and has a certain color palette and color breaks, values at which color changes
- the background of the map shows curved lines with constant latitude or longitude (graticule)
- the axis ticks show the latitude and longitude values
Polygons are a particular form of geometry; spatial geometries (points, lines, polygons, pixels) are discussed in detail in chapter 3. Polygons consist of sequences of points, connected by straight lines. How point locations of spatial data are expressed, or measured, is discussed in chapter 2. As can be seen from figure 1.1, lines of equal latitude and longitude do not form straight lines, indicating that some form of projection took place before plotting; projections are also discussed in chapter 2 and section 9.1.
The color values in figure 1.1 are derived
from numeric values of a variable,
BIR74, which has a
single value associated with each geometry or feature. Chapter
5 discusses such feature attributes, and the
way they can relate to feature geometries. In this case,
refers to birth counts, meaning counts over the region. This
implies that the count does not refer to a value associated with
every point inside the polygon, which the continuous color might
suggest, but rather measures an integral (sum) over the polygon.
# Simple feature collection with 100 features and 3 fields # Geometry type: MULTIPOLYGON # Dimension: XY # Bounding box: xmin: -84.3 ymin: 33.9 xmax: -75.5 ymax: 36.6 # Geodetic CRS: NAD27 # # A tibble: 100 × 4 # AREA BIR74 SID74 geom # <dbl> <dbl> <dbl> <MULTIPOLYGON [°]> # 1 0.114 1091 1 (((-81.5 36.2, -81.5 36.3, -81.6 36.3, -81.6 36.3, -81.7 36… # 2 0.061 487 0 (((-81.2 36.4, -81.2 36.4, -81.3 36.4, -81.3 36.4, -81.3 36… # 3 0.143 3188 5 (((-80.5 36.2, -80.5 36.3, -80.5 36.3, -80.5 36.3, -80.6 36… # # … with 97 more rows
The printed output shows:
- the (selected) dataset has 100 features (records) and 3 fields (attributes)
- the geometry type is
- it has dimension
XY, indicating that each point will consist of 2 coordinate values
- the range of x and y values of the geometry
- the coordinate reference system (CRS) is geodetic, with coordinates in degrees longitude and latitude associated to the
NAD27datum (chapter 2)
- the three selected attribute variables are followed by a variable
MULTIPOLYGONwith unit degrees that contains the polygon information
More complicated plots can involve facet plots with a map in each facet, as shown in figure 1.2.
year_labels = c("SID74" = "1974 - 1978", "SID79" = "1979 - 1984") nc.32119 %>% select(SID74, SID79) %>% pivot_longer(starts_with("SID")) -> nc_longer ggplot() + geom_sf(data = nc_longer, aes(fill = value)) + facet_wrap(~ name, ncol = 1, labeller = labeller(name = year_labels)) + scale_y_continuous(breaks = 34:36) + scale_fill_gradientn(colors = sf.colors(20)) + theme(panel.grid.major = element_line(color = "white"))
An interactive, leaflet-based map is obtained in figure ??.
In figure 1.1, the grey lines denote the graticule, a grid with lines along constant latitude or longitude. Clearly, these lines are not straight, which indicates that a projection of the data was used for which the x and y axis do not align with longitude and latitude. In figure ?? we see that the north boundary of North Carolina is plotted as a straight line again, indicating that another projection was used.
The ellipsoidal coordinates of the graticule of figure 1.1 are associated with a particular datum (here: NAD27), which implicates a set of rules what the shape of the Earth is and how it is attached to the Earth (to which point of the Earth is the origin associated, and how is it directed). If one would measure coordinates with a GPS device (e.g. a mobile phone) it would typically report coordinates associated with the WGS84 datum, which can be around 30 m different from the identical coordinate values when associated with NAD27.
Projections describe how we go back and forth between
ellipsoidal coordinates which are expressed as degrees latitude and longitude, pointing to locations on a shape approximating the Earth’s shape (an ellipsoid or spheroid), and
projected coordinates which are coordinates on a flat, two-dimensional coordinate system, used when plotting maps.
Datums transformations are associated with moving from one datum to another. Both topics are covered by spatial reference systems are described in more detail in chapter 2.
Polygon, point and line geometries are examples of vector data: point coordinates describe the “exact” locations that can be anywhere. Raster data on the other hand describe data where values are aligned on a raster, meaning on a regularly laid out lattice of usually square pixels. An example is shown in figure 1.3.
library(stars) par(mfrow = c(2, 2)) par(mar = rep(1, 4)) tif <- system.file("tif/L7_ETMs.tif", package = "stars") x <- read_stars(tif)[,,,1] image(x, main = "(a)") image(x[,1:10,1:10], text_values = TRUE, border = 'grey', main = "(b)") image(x, main = "(c)") set.seed(131) pts = st_sample(st_as_sfc(st_bbox(x)), 3) plot(pts, add = TRUE, pch = 3, col = 'blue') image(x, main = "(d)") plot(st_buffer(pts, 500), add = TRUE, pch = 3, border = 'blue', col = NA, lwd = 2)
Vector and raster data can be combined in different ways; for instance we can query the raster at the three points of figure 1.3(c),
# Simple feature collection with 3 features and 1 field # Geometry type: POINT # Dimension: XY # Bounding box: xmin: 290000 ymin: 9110000 xmax: 292000 ymax: 9120000 # Projected CRS: SIRGAS 2000 / UTM zone 25S # L7_ETMs.tif geometry # 1 80 POINT (290830 9114499) # 2 58 POINT (290019 9119219) # 3 63 POINT (291693 9116038)
or compute an aggregate, such as the average, over arbitrary regions such as the circles shown in figure 1.3(d):
# Simple feature collection with 3 features and 1 field # Geometry type: POLYGON # Dimension: XY # Bounding box: xmin: 290000 ymin: 9110000 xmax: 292000 ymax: 9120000 # Projected CRS: SIRGAS 2000 / UTM zone 25S # V1 geometry # 1 77.2 POLYGON ((291330 9114499, 2... # 2 60.1 POLYGON ((290519 9119219, 2... # 3 71.6 POLYGON ((292193 9116038, 2...
Other raster-to-vector conversions are discussed in 7.5 and include:
- converting raster pixels into point values
- converting raster pixels into small polygons, possibly merging polygons with identical values (“polygonize”)
- generating lines or polygons that delineate continuous pixel areas with a certain value range (“contour”)
Vector-to-raster conversions can be as simple as rasterizing polygons, as shown in figure 1.4. Other, more general vector-to-raster conversions that may involve statistical modelling include:
Raster dimensions describe how the rows and columns relate to spatial coordinates. Figure 1.5 shows a number of different possibilities.
x = 1:5 y = 1:4 d = st_dimensions(x = x, y = y, .raster = c("x", "y")) m = matrix(runif(20),5,4) r1 = st_as_stars(r = m, dimensions = d) r = attr(d, "raster") r$affine = c(0.2, -0.2) attr(d, "raster") = r r2 = st_as_stars(r = m, dimensions = d) r = attr(d, "raster") r$affine = c(0.1, -0.3) attr(d, "raster") = r r3 = st_as_stars(r = m, dimensions = d) x = c(1, 2, 3.5, 5, 6) y = c(1, 1.5, 3, 3.5) d = st_dimensions(x = x, y = y, .raster = c("x", "y")) r4 = st_as_stars(r = m, dimensions = d) grd = st_make_grid(cellsize = c(10,10), offset = c(-130,10), n = c(8,5), crs = st_crs(4326)) r5 = st_transform(grd, "+proj=laea +lon_0=-70 +lat_0=35") par(mfrow = c(2,3), mar = c(0.1, 1, 1.1, 1)) r1 = st_make_grid(cellsize = c(1,1), n = c(5,4), offset = c(0,0)) plot(r1, main = "regular") plot(st_geometry(st_as_sf(r2)), main = "rotated") plot(st_geometry(st_as_sf(r3)), main = "sheared") plot(st_geometry(st_as_sf(r4, as_points = FALSE)), main = "rectilinear") plot(st_geometry((r5)), main = "curvilinear")
Regular rasters like shown in figure 1.5 have a constant, not necessarily square cell size and axes aligned with the x and y (Easting and Northing) axes. Other raster types include those where the axes are no longer aligned with x and y (rotated), where axes are no longer perpendicular (sheared), or where cell size varies along a dimension (rectilinear). Finally, curvilinear rasters have cell size and/or direction properties that are no longer independent from the other raster dimension.
When a raster that is regular in a given coordinate reference system is projected to another raster while keeping each raster cell in tact, it changes shape and may become rectilinear (e.g. when going from ellipsoidal coordinates to Mercator, as in figure ??) or curvilinear (e.g. when going from ellipsoidal coordinates to Lambert Conic Conformal, as in figure 1.1). When reverting this procedure, one can recover the exact original raster.
Creating a new, regular grid in the new projection is called raster (or image) reprojection or warping (section 7.7). This process is lossy, irreversible, and may need to be informed whether raster cells should be interpolated, averaged or summed, whether they denote categorical variables, or whether resampling using nearest neighbours should be used; see also section 1.6.
A lot of spatial data is not just spatial, but in addition temporal. Just like any observation is associated with an observation location, it is associated with an observation time or period. The dataset on the North Carolina counties shown above contains disease cases counted over two time periods, shown in figure 1.2. Although the original dataset has these variables in two different columns, for plotting them these columns had to be stacked first, while repeating the associated geometries - a form called tidy by (Wickham 2014b). When we have longer time series associated with geometries, neither option - distributing time over multiple columns, or stacking columns while repeating geometries - works well, and a more effective way of storing such data would be a matrix or array, where one dimension refers to time, and the other(s) to space. The natural way for image or raster data is already to store them in matrices; time series of rasters then lead to a three-dimensional array. The general term for such data is a (spatiotemporal) data cube, where cube refers to arrays with any number of dimensions. Data cubes can refer to both raster and vector data, examples are given in chapter 6.
When we have spatial data with geometries that are not points but collections of points (multi-points, lines, polygons, pixels), then the attributes associated with these geometries has one of several different relationships to them. Attributes can have:
- a constant value for every point of the geometry
- a value that is unique to only this geometry, describing its identity
- a single value that is an aggregate over all points of the geometry
An example of a constant is land use or bedrock type of a polygon. An example of an identity is a county name. An example of an aggregate is the number of births over a given period of time, of a county.
The area with to which an attribute value refers to is called its support: aggregate properties have “block” (or polygon, or line) support, constant properties have “point” support (they apply to every point). Support matters when we manipulate the data. For instance, figure 1.4 was derived from a variable that has polygon support: the number of births per county. Rasterizing these values gives pixels with values that are associated to counties. The result of the rasterization is a meaningless map: the numeric values (“birth totals”) are not associated with the raster cells, and the county boundaries are no longer present. Totals of birth for the whole state can no longer be recovered from the pixel values. Ignoring support can easily lead to meaningless results. Chapter 5 discusses this further.
Raster cell values may have point support, e.g. when the cell records the elevation of the point at the cell centre in a digital elevation model, or cell support, e.g. when a satellite image pixel gives the color value averaged over (an area similar to the) pixel. Most file formats do not provide this information, yet it may be important to know when aggregating, regridding or warping rasters (section 7.7).
Although this book largely uses R and R packages for spatial data
science, a number of these packages use software libraries that were
not developed for R specifically. As an example, the dependency
of R package
sf on other R packages and system libraries is shown
in figure 1.6.
The C or C++ libraries used (GDAL, GEOS, PROJ, liblwgeom, s2geometry, NetCDF, udunits2) are all developed, maintained and used by (spatial) data science communities that are large and mostly different from the R community. By using these libraries, R users share how we understand what we are doing with these other communities. Because R, Python and Julia provide interactive interfaces to this software, many users get closer to these libraries than do users of other software based on these libraries. The first part of this book describes many of the concepts implemented in these libraries, which is relevant to spatial data science in general.
GDAL (“Geospatial Data Abstraction Library”) can be seen as the Swiss army knife of spatial data; besides for R it is being used in Python, QGIS, PostGIS, and more than 100 other software projects.
GDAL is a “library of libraries” – in order to read all these data sources it needs a large number of other libraries. It typically links to over 100 other libraries, each of which provides access to e.g. a particular data file format, database access, or web service.
Binary R packages distributed by CRAN contain only statically linked
code: CRAN does not want to make any assumptions about presence
of third-party libraries on the host system. As a consequence,
sf package is installed in binary form from CRAN, it
includes a copy of all the required external libraries as well as
their dependencies, which may amount to 100 Mb.
PROJ (or PR\(\phi\)J) is a library for cartographic projections and datum transformations: it converts spatial coordinates from one coordinate reference system to another. It comes with a large database of known projections and access to datum grids (high-precision pre-calculated values for datum transformations). It aligns with an international standard for coordinate reference systems (Lott 2015). Chapter 2 deals with coordinate systems, and PROJ.
GEOS (“Geometry Engine Open Source”) and s2geometry are two libraries for geometric operations. They are used to find measures (length, area, distance), and calculate predicates (do two geometries have any points in common?) or new geometries (which points do these two geometries have in common?). GEOS does this for flat, two-dimensional space (indicated by \(R^2\)), s2geometry does this for geometries on the sphere (indicated by \(S^2\)). Chapter 2 introduces coordinate reference systems, and chapter 4 discusses more about the differences between working with these two spaces.
NetCDF (UCAR 2020) refers to a file format as well as a C library for reading and writing NetCDF files. It allows the definition of arrays of any dimensionality, and is widely used for spatial and spatiotemporal information, especially in the (climate) modelling communities. Udunits2 (UCAR 2014; Pebesma et al. 2022) is a database and software library for units of measurement that allows the conversion of units, handles derived units, and supports user-defined units. The liblwgeom “library” is a software component of PostGIS (Obe and Hsu 2015) that contains several routines missing from GDAL or GEOS, including convenient access to GeographicLib routines (Karney 2013) that ship with PROJ.
- List five differences between raster and vector data.
- In addition to those listed below figure 1.1, list five further graphical components that are often found on a map.
- In your own words, why is the numeric information shown in figure 1.4 misleading (or meaningless)?
- Under which conditions would you expect strong differences when doing geometrical operations on \(S^2\), compared to doing them on \(R^2\)?
Karney, Charles FF. 2013. “Algorithms for Geodesics.” Journal of Geodesy 87 (1): 43–55. https://link.springer.com/content/pdf/10.1007/s00190-012-0578-z.pdf.
Lott, Roger. 2015. “Geographic Information-Well-Known Text Representation of Coordinate Reference Systems.” Open Geospatial Consortium. http://docs.opengeospatial.org/is/12-063r5/12-063r5.html.
Obe, Regina O, and Leo S Hsu. 2015. PostGIS in Action. Manning Publications Co.
Pebesma, Edzer, Thomas Mailund, Tomasz Kalinowski, and Iñaki Ucar. 2022. Units: Measurement Units for R Vectors. https://github.com/r-quantities/units/.
UCAR. 2014. UDUNITS 2.2.26 Manual. https://www.unidata.ucar.edu/software/udunits/udunits-current/doc/udunits/udunits2.html.
UCAR. 2020. The Netcdf User’s Guide. https://www.unidata.ucar.edu/software/netcdf/docs/user_guide.html.
Wickham, Hadley. 2014b. “Tidy Data.” Journal of Statistical Software 59 (1). https://www.jstatsoft.org/article/view/v059i10.
Wickham, Hadley, and Garret Grolemund. 2017. R for Data Science. O’Reilly. http://r4ds.had.co.nz/.