Converting point data

We’ve already encountered point data, so we’ll start with that.

We can do this with our city data from earlier in the semester. The key function here is st_as_sf, which takes an ordinary dataframe and turns it into an sf object. In order for that to happen, you need to tell it which columns are geographic, so that it can bind those into a new geometry feature; and a “crs” (coordinate reference system) which specifies what kind of coordinates you’re passing in. The CRS argument is required, but will almost always be 4326, which is the EPSG code to represent plain latitude and longitude. EPSG stands for “European Petroleum Survey Group;” it is a no-longer-existing (merged into the International Association of Oil & Gas Producers) industrial group that wanted to keep lists of different map projections straight.

The code below creates a spatial features frame from a set of cities

If you type cities in, you’ll see a variety of different metadata, including that this is a point set: that it has a bbox (a bounding box); and two pieces of information about the map projection. One is the “EPSG” code; the other is a ‘proj4’ string that’s another way of setting the projection. We’ll talk more about projections later in this chapter.

The rest is simply the same as our dataset, except that it now has one additional column called ‘geometry.’

      cities <- CESTA |>
  # Removes both NA values and those coded as 0.
  filter(!is.na(LAT), LAT != 0) |>
  st_as_sf(coords = c("LON", "LAT"), crs = 4326)

cities
    

We can plot this frame by using a new geom called geom_sf. Note that we no longer need to include ‘x’ and ‘y’: the coordinate system takes care of that for us.

      cities |> filter(ST != "AK", ST!="HI") |>
ggplot() + geom_sf()
    

Any of the other aesthetics we know from ggplot2 can be applied. size changes the size of points; fill changes their colors; and the use of geom_sf_label makes it possible to add a set of labels on top of a chart.

Note an additional ggplot trick here for plotting multiple datasets at the same time. Inside geom_sf_label, we use a different dataset: cities |> filter(1990> 1000000). Without that parameter, we would be trying to label all the cities in the dataset.

With mapping in R, it is frequently very useful to plot multiple items on top of each other. You’ll see this same pattern again and again in this notebook.

      cities |> 
  filter(ST != "AK", ST!="HI", `1980` > 0, `1990` > 0) |>
  ggplot() + 
  geom_sf(alpha = 0.25) +
  geom_sf_label(data = cities |> filter(`1990` > 1000000)) +
  aes(size = `1990`, color=ST, label = City) + 
  theme(legend.position="none") +
  scale_size_continuous(trans='sqrt') + 
  labs(title="Eight US cities had populations above 1,000,000 in 1990.")
    

Reading shapefiles.

There are several different ways to make an sf object. One is to just read in any shapefile, geojson, or other spatial object using the read_sf function. This is essentially the same as the read_csv functions with which you’re already familiar.

(One niggling note: you might find old code online that uses st_read instead of read_sf. This is less likely to work with tidyverse packages, and so you should always use read_sf.)

      us_states_file_location = system.file("extdata/states/cb_2018_us_state_20m.shp", package="HumanitiesDataAnalysis")
states = read_sf(us_states_file_location)
ggplot(states) + geom_sf()
    

Choosing projections.

If you are plotting either the whole world, or some relatively obscure part of it, though, you may need to make a decision. Since the Earth is (almost) a sphere, any flat projection will bring some kind of distortion. The choice of map projection is about meeting two incompatible goals. When you wade through options, consider the tradeoffs you make on the following fronts:

1. Preserving area. A projection that preserves area is called an ‘equal area’ projection. A good choice at the global scale is the ‘Mollweide’ projection; for national mapping of large countries (like Russia, the US, or China) an Albers projection is best. At the local scale, such as a single city, area distortions are unlikely to be important and you can focus on the other goal:

2. Preserving shape. A projection that preserves local shapes is called a ‘conformal’ projection. At the global scale, the Mercator projection is widely used. (Too many well-educated people know only one thing about about map projections; that the Mercator projection is distortionary and therefore bad.) For regional mapping, the Lambert conformal cylindrical projection is frequently deployed in regions far from the equator; closer to the equator, the Mercator projection is fine. For geographic entities that extend mostly north-south (like Chile or California) a transverse Mercator projection works better than the traditional Mercator.

We need multiple ways to store information.

This explains why we need to choose a coordinate system for plotting. But one of the things that most perplexes newcomers to digital mapping is that

In a minute, we will talk about spatial joins. But the most basic join between our two sets fails, if you try to run it, because they are in different coordinate systems before being plotted.

      cities |> st_join(states)
    

This causes no end of trouble for beginners to mapping. Why not just store everything as longitude and latitude, and only worry about projections when it’s time to actually draw a map? Coordinate systems are a pain because the earth is not an even cartesian plane. There is no way around this. And it matters not just for plotting, but for storing data. There are three reasons.

To think about this, imagine that we are drawing not all of America but just–say–a map of Washington Square Park in New York City, including pathways, underground sewage, benches, etc. Longitude and latitude would fail for three reasons:

1. The numbers would get absurdly small; a bench might be 0.0005 degrees of longitude long.
2. The earth is not a sphere; in fact, it is a slightly lumpy ball of rock. A circle drawn around the equator is longer than one drawn around the two poles; and different parts of the earth are slightly higher or lower than others. The overall shape of the earth is known as the ‘geoid’; but since this is an ongoign project of measurement, most projections assume a somewhat simplified shape that will be more appropriate in some parts of the world than others.
3. The continents themselves are not fixed in location. Australia is moving at 7cm a year; if you were to precisely map the locations of a park in Sydney, in twenty years all the locations would be off by over a meter. That might not matter for a web map, but it could be disastrous for an archeological dig.

For all these reasons, most local data uses a local projection city. New York City, for example, generally uses a version of the New York State Plane. The same holds for most other regions. There are two important ways to define a projection. One is as a number; this refers to a list of “EPSG” codes. You’ll generally come across them if you search the Internet. If you search, say, for “EPSG Colorado map projection,” it will lead you to the page for EPSG:26954; pasting this number in will then give you a good projection for Colorado that captures how the state is narrower at the top than the bottom even though i’s ostensibly a square.

      cities |>
  filter(ST=="CO") |>
  st_transform(crs = 26954) |>
  ggplot() +
  geom_sf(data = states |> filter(STUSPS == "CO") |> st_transform(crs=26954)) + 
  geom_sf() + theme_minimal()# + geom_sf_label(aes(label = City))
    

Data manipulation

You can treat sf objects just like dataframes for most purposes; for instance, you can mutate countries to create a field for per-capita GDP and then make a boxplot to understand the categories it uses. This is not a map, though it’s made from map data. This is worth remembering–too often, people insist of visualizing geographic data as maps when another form would do better.

      countries |>
  mutate(gdp_per_capita = gdp_md_est / pop_est) |>
  ggplot() + aes(x = str_wrap(economy, 15), y = gdp_per_capita) + geom_boxplot() + scale_y_log10() +
  geom_point(aes(size = pop_est, color = continent), position = position_jitter()) +
  theme(legend.position = "bottom")
    

Cartography

But you can also view sf features as a map. The easiest way to do this is by using the geom_sf function in ggplot. Just as geom_point adds a point, geom_sf adds any spatial element to a chart.

      countries |>
 # st_transform(crs = "+proj=moll") |>
  ggplot() +
  geom_sf(aes(fill = economy), lwd=0.1) + 
  theme_minimal() +
  labs(title="World Countries by economic classification")
    

To change the color of an object, use st_transform with the crs that you want. A Mollweide projection can be created using a proj4string that simply says moll; it’s often a decent choice for mapping the globe if you wish to preserve equal area.

      #make_background = function(crs) {st_graticule(ndiscr = 1000, margin = 10e-5) |>
#  st_transform(crs = crs) |>
#  st_convex_hull() |>
#  summarise(geometry = st_union(geometry))
#}
rnaturalearthdata::sovereignty50 |>
  st_as_sf() |>
  st_transform(crs = "+proj=moll") |>
  ggplot() +
  geom_sf(aes(fill = economy), lwd=.15)
    

Spatial Summaries

Most of the tidyverse functions we use will continue to work on spatial objects. For instance, a summary can still be used to get the counts by time period. But summary has an additional action that it performs on geometries.

Here’s an example. If we extract villas from the Pleiades dataset, we haveover 2000 points; but if we summarize by time period, we can take just the first 8 items. But if we plot it, we can see that each of these 8 features has dozens or hundreds of individual points inside of it. This can be extraordinarily powerful.

      villas = pleiades |> collect() |>
 unnest(featureTypes) |> filter(featureTypes == "villa") 
villas |> select(featureTypes)
    

      villa_counts = villas |>
  unnest(timePeriodsKeys) |>
  group_by(timePeriodsKeys) |>
  summarize(count = n()) |>
  arrange(desc(count)) |> head(8)

villa_counts |> ggplot() + geom_sf(data=provinces) + geom_sf() + facet_wrap(~timePeriodsKeys) 
    

By using additional transformations, this can create different types of features. For example, the function st_convex_hull draws a line around all the entries in a point: if we do that on our collections of cities, we get blobs that approximate the shapes of the states. This can be useful if you’re working with data that doesn’t have a defined geometry.

      cities |> group_by(ST) |> 
  summarize(count = n(), urbanPop2010 = sum(`2010`)) |> 
  st_transform(crs = 2163) |>
  st_convex_hull() |> ggplot() + geom_sf() + 
  aes(fill = urbanPop2010/count) + 
  scale_fill_viridis_c(trans="log", labels = scales::comma)
    

Spatial joins

In some historical cases, the outlines may be all that you have. But in general, spatial joins allow a better way to plot areas than janky estimates like the one above.

Suppose we want to plot the number of villas in each province. This requires locating each point inside a polygon. A spatial join takes two sf features and links them when they share the same space.

By default, two items are linked if they intersect; but you can also use a variety of other functions, including requiring that all the points in one object be within those on the right, that they contain, that they be within a certain distance, and so forth. (See ?st_join for a full list.) If, for instance, you wanted to find all villas within 10 miles of a Roman road, you would use villas |> st_join(roads, st_is_within_distance, unit).

Two poihts are important to keep in mind when doing a join:

1. The geometries returned will be from the left side of the join. So if you have a set of points from Pleiades and set of provinces, a join which starts with the points will return more points; a join that starts with the provinces will return more provinces.

2. An st_join operation can create many duplicate points, so it is best to do any counting operations before performing the join. For the most part in this text we ignore these issues of computational complexity; but shapefiles can be quite elaborate. Sometimes this means that you’ll actually want to do two joins.

      villas = pleiades |> 
  unnest(featureTypes) |>
  filter(featureTypes |> str_detect("villa")) |>
  st_join(provinces) |>
  count(name, name = "Number of villas")

provinces |> st_join(villas) |>
  ggplot() + geom_sf(aes(fill=`Number of villas`)) + 
  scale_fill_viridis_c("Number of villas in province")
    

There are also some functions that can be used inside a mutate call. One is st_area, which calculates the area of a region. For villas, it probably makes more sense to plot in terms of ‘villas per square kilometer’; the following code shows that happening.

      provinces |>
  st_join(villas) |>
  mutate(area = st_area(geometry)) |>
  # as.numeric(area) is necessary to avoid an error
  # since area is cast as 'units', not number.
  ggplot() + geom_sf() + aes(fill = `Number of villas` / as.numeric(area) * 1e06) + 
  scale_fill_viridis_c("Number of villas per square kilometer")
    

Grid plotting

One last, advanced geo-computation element worth considering is using not previously existing shapefiles, but instead creating your own using a grid. If you want, for example, to see how many total points there are, you can bin them into a spatial grid before plotting, by first creating a grid using st_make_grid on a bounding box derived from your data.

      bbox <- countries |> filter(admin %in%
  c("Ireland", "Iran", "Egypt", "Morocco", "Estonia")) |>
  EUproject() |> st_bbox()

grid <- st_make_grid(bbox |> st_as_sfc(), n = c(60, 60), square = FALSE) |> st_sf() |> mutate(grid_id = 1:n())

counts = pleiades |>
  st_join(grid) |>
  as_tibble() |>
  select(-geometry) |>
  group_by(grid_id) |>
  summarize(points = n())

grid |>
  inner_join(counts) |>
  filter(points > 20) |>
  ggplot() +
  scale_fill_viridis_c(trans = "log") +
  geom_sf(data = ancient |> summarize(shape = 1), fill = "brown", color = "#00000000") +
  theme(panel.background = element_rect(fill = "aliceblue")) + geom_sf(aes(fill = points)) + coord_sf(expand = FALSE) + 
  labs(title = "Count")
    

Looking at the pure counts tells us something. I would not have known ahead of time that the area around Carthage is as heavily represented as around Rome or Athens. But we can use these as documents instead, as well.

      types |>
  unnest(featureTypes) |>
  st_join(grid) |>
  count(grid_id, featureTypes) |>
  arrange(-n) |>
  group_by(grid_id) |>
  filter(!featureTypes %in% c("labeled feature", " settlement-modern", " unlabeled", "unknown")) -> counts

tops = counts |>
  summarize_llr(featureTypes, count = n) |>
  group_by(grid_id) |>
  filter(.dunning_llr > .05) |>
  arrange(-.dunning_llr) |>
  slice(1) |>
  ungroup() |>
  arrange(-.dunning_llr) |>
  ungroup() |>
  mutate(feature = fct_lump(featureTypes, 7)) |>
  filter(!feature %in% c("Other", "people"))

st_join(grid, tops) |>
  ggplot() + scale_fill_brewer(type = "qual") +
  geom_sf(data = ancient |> summarize(shape = 1), fill = "brown", alpha = .1, color = "#00000000") +
  theme(panel.background = element_rect(fill = "aliceblue")) +
  geom_sf(aes(fill = feature, alpha = .dunning_llr), color = "#00000000") + coord_sf(expand = FALSE) + scale_alpha_continuous(range = c(.5, 1))
    

      rnaturalearthdata::countries110 |>
  st_as_sf() |>
  ggplot() + geom_sf(aes(fill=type))