3.2 Vector attribute manipulation

Geographic vector data in R are well supported by sf, a class which extends the data.frame.Thus sf objects have one column per attribute variable (such as ‘name’) and one row per observation, or feature (e.g., per bus station).sf objects also have a special column to contain geometry data, usually named geometry.The geometry column is special because it is a list column, which can contain multiple geographic entities (points, lines, polygons) per row.This was described in Chapter 2, which demonstrated how generic methods such as plot() and summary() work on sf objects.sf also provides methods that allow sf objects to behave like regular data frames, as illustrated by other sf-specific methods that were originally developed for data frames:

  1. methods(class = "sf") # methods for sf objects, first 12 shown
  1. #>  [1] aggregate             cbind                 coerce               
  2. #>  [4] initialize            merge                 plot                 
  3. #>  [7] print                 rbind                 [                    
  4. #> [10] [[<-                  $<-                   show

Many of these functions, including rbind() (for binding rows of data together) and $<- (for creating new columns) were developed for data frames.A key feature of sf objects is that they store spatial and non-spatial data in the same way, as columns in a data.frame.

The geometry column of sf objects is typically called geometry but any name can be used.The following command, for example, creates a geometry column named g:

st_sf(data.frame(n = world$name_long), g = world$geom)
This enables geometries imported from spatial databases to have a variety of names such as wkb_geometry and the_geom.

sf objects also support tibble and tbl classes used in the tidyverse, allowing ‘tidy’ data analysis workflows for spatial data.Thus sf enables the full power of R’s data analysis capabilities to be unleashed on geographic data.Before using these capabilities it is worth re-capping how to discover the basic properties of vector data objects.Let’s start by using base R functions to get a measure of the world dataset:

  1. dim(world) # it is a 2 dimensional object, with rows and columns
  2. #> [1] 177  11
  3. nrow(world) # how many rows?
  4. #> [1] 177
  5. ncol(world) # how many columns?
  6. #> [1] 11

Our dataset contains ten non-geographic columns (and one geometry list column) with almost 200 rows representing the world’s countries.Extracting the attribute data of an sf object is the same as removing its geometry:

  1. world_df = st_drop_geometry(world)
  2. class(world_df)
  3. #> [1] "data.frame"

This can be useful if the geometry column causes problems, e.g., by occupying large amounts of RAM, or to focus the attention on the attribute data.For most cases, however, there is no harm in keeping the geometry column because non-spatial data operations on sf objects only change an object’s geometry when appropriate (e.g., by dissolving borders between adjacent polygons following aggregation).This means that proficiency with attribute data in sf objects equates to proficiency with data frames in R.

For many applications, the tidyverse package dplyr offers the most effective and intuitive approach for working with data frames.Tidyverse compatibility is an advantage of sf over its predecessor sp, but there are some pitfalls to avoid (see the supplementary tidyverse-pitfalls vignette at geocompr.github.io for details).

3.2.1 Vector attribute subsetting

Base R subsetting functions include [, subset() and $.dplyr subsetting functions include select(), filter(), and pull().Both sets of functions preserve the spatial components of attribute data in sf objects.

The [ operator can subset both rows and columns.You use indices to specify the elements you wish to extract from an object, e.g., object[i, j], with i and j typically being numbers or logical vectors — TRUEs and FALSEs — representing rows and columns (they can also be character strings, indicating row or column names).Leaving i or j empty returns all rows or columns, so world[1:5, ] returns the first five rows and all columns.The examples below demonstrate subsetting with base R.The results are not shown; check the results on your own computer:

  1. world[1:6, ] # subset rows by position
  2. world[, 1:3] # subset columns by position
  3. world[, c("name_long", "lifeExp")] # subset columns by name

A demonstration of the utility of using logical vectors for subsetting is shown in the code chunk below.This creates a new object, small_countries, containing nations whose surface area is smaller than 10,000 km2:

  1. sel_area = world$area_km2 < 10000
  2. summary(sel_area) # a logical vector
  3. #>    Mode   FALSE    TRUE 
  4. #> logical     170       7
  5. small_countries = world[sel_area, ]

The intermediary sel_area is a logical vector that shows that only seven countries match the query.A more concise command, which omits the intermediary object, generates the same result:

  1. small_countries = world[world$area_km2 < 10000, ]

The base R function subset() provides yet another way to achieve the same result:

  1. small_countries = subset(world, area_km2 < 10000)

Base R functions are mature and widely used.However, the more recent dplyr approach has several advantages.It enables intuitive workflows.It is fast, due to its C++ backend.This is especially useful when working with big data as well as dplyr’s database integration.The main dplyr subsetting functions are select(), slice(), filter() and pull().

raster and dplyr packages have a function called select(). When using both packages, the function in the most recently attached package will be used, ‘masking’ the incumbent function. This can generate error messages containing text like: unable to find an inherited method for function ‘select’ for signature ‘“sf”’. To avoid this error message, and prevent ambiguity, we use the long-form function name, prefixed by the package name and two colons (usually omitted from R scripts for concise code): dplyr::select().

select() selects columns by name or position.For example, you could select only two columns, name_long and pop, with the following command (note the sticky geom column remains):

  1. world1 = dplyr::select(world, name_long, pop)
  2. names(world1)
  3. #> [1] "name_long" "pop"       "geom"

select() also allows subsetting of a range of columns with the help of the : operator:

  1. # all columns between name_long and pop (inclusive)
  2. world2 = dplyr::select(world, name_long:pop)

Omit specific columns with the - operator:

  1. # all columns except subregion and area_km2 (inclusive)
  2. world3 = dplyr::select(world, -subregion, -area_km2)

Conveniently, select() lets you subset and rename columns at the same time, for example:

  1. world4 = dplyr::select(world, name_long, population = pop)
  2. names(world4)
  3. #> [1] "name_long"  "population" "geom"

This is more concise than the base R equivalent:

  1. world5 = world[, c("name_long", "pop")] # subset columns by name
  2. names(world5)[names(world5) == "pop"] = "population" # rename column manually

select() also works with ‘helper functions’ for advanced subsetting operations, including contains(), starts_with() and num_range() (see the help page with ?select for details).

Most dplyr verbs return a data frame.To extract a single vector, one has to explicitly use the pull() command.The subsetting operator in base R (see ?[), by contrast, tries to return objects in the lowest possible dimension.This means selecting a single column returns a vector in base R.To turn off this behavior, set the drop argument to FALSE.

  1. # create throw-away data frame
  2. d = data.frame(pop = 1:10, area = 1:10)
  3. # return data frame object when selecting a single column
  4. d[, "pop", drop = FALSE] # equivalent to d["pop"]
  5. select(d, pop)
  6. # return a vector when selecting a single column
  7. d[, "pop"]
  8. pull(d, pop)

Due to the sticky geometry column, selecting a single attribute from an sf-object with the help of [() returns also a data frame.Contrastingly, pull() and $ will give back a vector.

  1. # data frame object
  2. world[, "pop"]
  3. # vector objects
  4. world$pop
  5. pull(world, pop)

slice() is the row-equivalent of select().The following code chunk, for example, selects the 3rd to 5th rows:

  1. slice(world, 3:5)

filter() is dplyr’s equivalent of base R’s subset() function.It keeps only rows matching given criteria, e.g., only countries with a very high average of life expectancy:

  1. # Countries with a life expectancy longer than 82 years
  2. world6 = filter(world, lifeExp > 82)

The standard set of comparison operators can be used in the filter() function, as illustrated in Table 3.1:

Table 3.1: Comparison operators that return Booleans (TRUE/FALSE).
SymbolName
==Equal to
!=Not equal to
>, <Greater/Less than
>=, <=Greater/Less than or equal
&, |, !Logical operators: And, Or, Not

dplyr works well with the ‘pipe’ operator %>%, which takes its name from the Unix pipe | (Grolemund and Wickham 2016).It enables expressive code: the output of a previous function becomes the first argument of the next function, enabling chaining.This is illustrated below, in which only countries from Asia are filtered from the world dataset, next the object is subset by columns (name_long and continent) and the first five rows (result not shown).

  1. world7 = world %>%
  2.   filter(continent == "Asia") %>%
  3.   dplyr::select(name_long, continent) %>%
  4.   slice(1:5)

The above chunk shows how the pipe operator allows commands to be written in a clear order:the above run from top to bottom (line-by-line) and left to right.The alternative to %>% is nested function calls, which is harder to read:

  1. world8 = slice(
  2.   dplyr::select(
  3.     filter(world, continent == "Asia"),
  4.     name_long, continent),
  5.   1:5)

3.2.2 Vector attribute aggregation

Aggregation operations summarize datasets by a ‘grouping variable’, typically an attribute column (spatial aggregation is covered in the next chapter).An example of attribute aggregation is calculating the number of people per continent based on country-level data (one row per country).The world dataset contains the necessary ingredients: the columns pop and continent, the population and the grouping variable, respectively.The aim is to find the sum() of country populations for each continent.This can be done with the base R function aggregate() as follows:

  1. world_agg1 = aggregate(pop ~ continent, FUN = sum, data = world, na.rm = TRUE)
  2. class(world_agg1)
  3. #> [1] "data.frame"

The result is a non-spatial data frame with six rows, one per continent, and two columns reporting the name and population of each continent (see Table 3.2 with results for the top 3 most populous continents).

aggregate() is a generic function which means that it behaves differently depending on its inputs.sf provides a function that can be called directly with sf:::aggregate() that is activated when a by argument is provided, rather than using the ~ to refer to the grouping variable:

  1. world_agg2 = aggregate(world["pop"], by = list(world$continent),
  2.                        FUN = sum, na.rm = TRUE)
  3. class(world_agg2)
  4. #> [1] "sf"         "data.frame"

As illustrated above, an object of class sf is returned this time.world_agg2 which is a spatial object containing 6 polygons representing the columns of the world.

summarize() is the dplyr equivalent of aggregate().It usually follows group_by(), which specifies the grouping variable, as illustrated below:

  1. world_agg3 = world %>%
  2.   group_by(continent) %>%
  3.   summarize(pop = sum(pop, na.rm = TRUE))

This approach is flexible and gives control over the new column names.This is illustrated below: the command calculates the Earth’s population (~7 billion) and number of countries (result not shown):

  1. world %>% 
  2.   summarize(pop = sum(pop, na.rm = TRUE), n = n())

In the previous code chunk pop and n are column names in the result.sum() and n() were the aggregating functions.The result is an sf object with a single row representing the world (this works thanks to the geometric operation ‘union’, as explained in Section 5.2.6).

Let’s combine what we have learned so far about dplyr by chaining together functions to find the world’s 3 most populous continents (with dplyr::top_n() ) and the number of countries they contain (the result of this command is presented in Table 3.2):

  1. world %>% 
  2.   dplyr::select(pop, continent) %>% 
  3.   group_by(continent) %>% 
  4.   summarize(pop = sum(pop, na.rm = TRUE), n_countries = n()) %>% 
  5.   top_n(n = 3, wt = pop) %>%
  6.   st_drop_geometry()
Table 3.2: The top 3 most populous continents, and the number of countries in each.
continentpopn_countries
Africa115494663351
Asia431140805947
Europe66903625639

More details are provided in the help pages (which can be accessed via ?summarize and vignette(package = "dplyr") and Chapter 5 of R for Data Science.

3.2.3 Vector attribute joining

Combining data from different sources is a common task in data preparation.Joins do this by combining tables based on a shared ‘key’ variable.dplyr has multiple join functions including left_join() and inner_join() — see vignette("two-table") for a full list.These function names follow conventions used in the database language SQL(Grolemund and Wickham 2016, Chapter 13); using them to join non-spatial datasets to sf objects is the focus of this section.dplyr join functions work the same on data frames and sf objects, the only important difference being the geometry list column.The result of data joins can be either an sf or data.frame object.The most common type of attribute join on spatial data takes an sf object as the first argument and adds columns to it from a data.frame specified as the second argument.

To demonstrate joins, we will combine data on coffee production with the world dataset.The coffee data is in a data frame called coffee_data from the spData package (see ?coffee_data for details).It has 3 columns:name_long names major coffee-producing nations and coffee_production_2016 and coffee_production_2017 contain estimated values for coffee production in units of 60-kg bags in each year.A ‘left join’, which preserves the first dataset, merges world with coffee_data:

  1. world_coffee = left_join(world, coffee_data)
  2. #> Joining, by = "name_long"
  3. class(world_coffee)
  4. #> [1] "sf"         "data.frame"

Because the input datasets share a ‘key variable’ (name_long) the join worked without using the by argument (see ?left_join for details).The result is an sf object identical to the original world object but with two new variables (with column indices 11 and 12) on coffee production.This can be plotted as a map, as illustrated in Figure 3.1, generated with the plot() function below:

  1. names(world_coffee)
  2. #>  [1] "iso_a2"                 "name_long"             
  3. #>  [3] "continent"              "region_un"             
  4. #>  [5] "subregion"              "type"                  
  5. #>  [7] "area_km2"               "pop"                   
  6. #>  [9] "lifeExp"                "gdpPercap"             
  7. #> [11] "coffee_production_2016" "coffee_production_2017"
  8. #> [13] "geom"
  9. plot(world_coffee["coffee_production_2017"])

World coffee production (thousand 60-kg bags) by country, 2017. Source: International Coffee Organization.
Figure 3.1: World coffee production (thousand 60-kg bags) by country, 2017. Source: International Coffee Organization.

For joining to work, a ‘key variable’ must be supplied in both datasets.By default dplyr uses all variables with matching names.In this case, both world_coffee and world objects contained a variable called name_long, explaining the message Joining, by = "name_long".In the majority of cases where variable names are not the same, you have two options:

  • Rename the key variable in one of the objects so they match.
  • Use the by argument to specify the joining variables.
    The latter approach is demonstrated below on a renamed version of coffee_data:
  1. coffee_renamed = rename(coffee_data, nm = name_long)
  2. world_coffee2 = left_join(world, coffee_renamed, by = c(name_long = "nm"))

Note that the name in the original object is kept, meaning that world_coffee and the new object world_coffee2 are identical.Another feature of the result is that it has the same number of rows as the original dataset.Although there are only 47 rows of data in coffee_data, all 177 the country records are kept intact in world_coffee and world_coffee2:rows in the original dataset with no match are assigned NA values for the new coffee production variables.What if we only want to keep countries that have a match in the key variable?In that case an inner join can be used:

  1. world_coffee_inner = inner_join(world, coffee_data)
  2. #> Joining, by = "name_long"
  3. nrow(world_coffee_inner)
  4. #> [1] 45

Note that the result of inner_join() has only 45 rows compared with 47 in coffee_data.What happened to the remaining rows?We can identify the rows that did not match using the setdiff() function as follows:

  1. setdiff(coffee_data$name_long, world$name_long)
  2. #> [1] "Congo, Dem. Rep. of" "Others"

The result shows that Others accounts for one row not present in the world dataset and that the name of the Democratic Republic of the Congo accounts for the other:it has been abbreviated, causing the join to miss it.The following command uses a string matching (regex) function from the stringr package to confirm what Congo, Dem. Rep. of should be:

  1. str_subset(world$name_long, "Dem*.+Congo")
  2. #> [1] "Democratic Republic of the Congo"

To fix this issue, we will create a new version of coffee_data and update the name.inner_join()ing the updated data frame returns a result with all 46 coffee-producing nations:

  1. coffee_data$name_long[grepl("Congo,", coffee_data$name_long)] = 
  2.   str_subset(world$name_long, "Dem*.+Congo")
  3. world_coffee_match = inner_join(world, coffee_data)
  4. #> Joining, by = "name_long"
  5. nrow(world_coffee_match)
  6. #> [1] 46

It is also possible to join in the other direction: starting with a non-spatial dataset and adding variables from a simple features object.This is demonstrated below, which starts with the coffeedata object and adds variables from the original world dataset.In contrast with the previous joins, the result is _not another simple feature object, but a data frame in the form of a tidyverse tibble:the output of a join tends to match its first argument:

  1. coffee_world = left_join(coffee_data, world)
  2. #> Joining, by = "name_long"
  3. class(coffee_world)
  4. #> [1] "tbl_df"     "tbl"        "data.frame"

In most cases, the geometry column is only useful in an sf object.The geometry column can only be used for creating maps and spatial operations if R ‘knows’ it is a spatial object, defined by a spatial package such as sf.Fortunately, non-spatial data frames with a geometry list column (like coffee_world) can be coerced into an sf object as follows: st_as_sf(coffee_world).

This section covers the majority joining use cases.For more information, we recommend Grolemund and Wickham (2016), the join vignette in the geocompkg package that accompanies this book, and documentation of the data.table package.15Another type of join is a spatial join, covered in the next chapter (Section 4.2.3).

3.2.4 Creating attributes and removing spatial information

Often, we would like to create a new column based on already existing columns.For example, we want to calculate population density for each country.For this we need to divide a population column, here pop, by an area column, here area_km2 with unit area in square kilometers.Using base R, we can type:

  1. world_new = world # do not overwrite our original data
  2. world_new$pop_dens = world_new$pop / world_new$area_km2

Alternatively, we can use one of dplyr functions - mutate() or transmute().mutate() adds new columns at the penultimate position in the sf object (the last one is reserved for the geometry):

  1. world %>% 
  2.   mutate(pop_dens = pop / area_km2)

The difference between mutate() and transmute() is that the latter skips all other existing columns (except for the sticky geometry column):

  1. world %>% 
  2.   transmute(pop_dens = pop / area_km2)

unite() pastes together existing columns.For example, we want to combine the continent and region_un columns into a new column named con_reg.Additionally, we can define a separator (here: a colon :) which defines how the values of the input columns should be joined, and if the original columns should be removed (here: TRUE):

  1. world_unite = world %>%
  2.   unite("con_reg", continent:region_un, sep = ":", remove = TRUE)

The separate() function does the opposite of unite(): it splits one column into multiple columns using either a regular expression or character positions.

  1. world_separate = world_unite %>% 
  2.   separate(con_reg, c("continent", "region_un"), sep = ":")

The two functions rename() and set_names() are useful for renaming columns.The first replaces an old name with a new one.The following command, for example, renames the lengthy name_long column to simply name:

  1. world %>% 
  2.   rename(name = name_long)

set_names() changes all column names at once, and requires a character vector with a name matching each column.This is illustrated below, which outputs the same world object, but with very short names:

  1. new_names = c("i", "n", "c", "r", "s", "t", "a", "p", "l", "gP", "geom")
  2. world %>% 
  3.   set_names(new_names)

It is important to note that attribute data operations preserve the geometry of the simple features.As mentioned at the outset of the chapter, it can be useful to remove the geometry.To do this, you have to explicitly remove it because sf explicitly makes the geometry column sticky.This behavior ensures that data frame operations do not accidentally remove the geometry column.Hence, an approach such as select(world, -geom) will be unsuccessful and you should instead use st_drop_geometry().16

  1. world_data = world %>% st_drop_geometry()
  2. class(world_data)
  3. #> [1] "data.frame"