06-Attributes.Rmd

# Feature attributes {#featureattributes}

Feature _attributes_ refer to the properties of features ("things")
that do not describe the feature's geometry. Feature attributes can
be _derived_ from geometry (e.g. length of a `LINESTRING`, area
of a `POLYGON`) but they can also refer to completely different
properties, such as 

* the name of a street or a county, 
* the number of people living in a country, 
* the type of a road
* the soil type in a polygon from a soil map.
* the opening hours of a shop
* the body weight of an animal
* the NO$_2$ concentration measured at an air quality monitoring station

Although temporal properties of features are no less fundamental
than their spatial properties, the simple feature access standard
and consequently the `sf` package does not give time a similar role
as space; more on that in chapter \@ref(raster).

Most `sf` objects will contain both geometries and attributes for
features. The geometric operations described in the previous chapter
(\@ref(geommanip)) operate on geometries _only_, and may occasionally
add attributes, but will not modify attributes present. 

In all these cases, attribute _values_ remain unmodified. At first
sight, that looks rather harmless. But if we look into a simple
case of replacing a county boundary with a county centroid, as in
```{r}
library(sf)
library(dplyr)
system.file("gpkg/nc.gpkg", package="sf") %>%
	read_sf() %>%
	st_transform(32119) %>%
	select(BIR74, SID74, NAME) %>%
	st_centroid() %>%
	head(n = 1) -> x # save as x
st_geometry(x)[[1]]
```
we receive a warning. This warning is justified for the first two
variables shown (total births and number of SID disease cases,
1974) which, as such, are _not associated with_ a feature whose
geometry is `POINT (385605.4 300303.5)`. The third variable,
`NAME` is however still the county name for the point indicated,
but the point geometry no longer _is_ the county geometry.

## Attribute-geometry relationships {#agr}

Changing the feature geometry without changing the feature attributes
does change the _feature_, since the feature is characterised by
the combination of geometry and attributes. Can we, ahead of time,
predict whether the resulting feature will still meaningfully relate
to the attribute data when we replace all geometries for instance
with their convex hull or centroid? It depends.

Take the example of a road, represented by a `LINESTRING`, which has
an attribute property _road width_ equal to 10 m. What can we say about
the road width of an arbitray subsectin of this road? That depends
on whether the attribute road length describes, for instance the
road width everywhere, meaning that road width is constant along the
road, or whether it describes an aggregate property, such as minimum
or average road width.  In case of the minimum, for an arbitrary
subsection of the road one could still argue that the minimum
road width must be at least as large as the minimum road width for
the whole segment, but it may no longer be _the minimum_ for that
subsection. This gives us two "types" for the attribute-geometry
relationship (AGR):

* **constant** the attribute value is valid everywhere in or over the geometry
* **aggregate** the attribute is an aggregate, a summary value over the geometry

For polygon data, typical examples of **constant** AGR are

* land use for a land use polygon
* rock units or geologic strata in a geological map
* soil type in a soil map
* elevation class in a elevation map that shows elevation as classes
* climate zone in a climate zone map

Typical examples for the **aggregate** AGR are

* population, either as number of persons or as population density
* other socio-economic variables, summarised by area
* total emission of pollutants by region
* block mean NO$_2$ concentrations, as e.g. obtained by block kriging or a dispersion model that predicts areal means

A third type of AGR is that where an attribute identifies a feature
geometry. The example above is county `NAME`: the name identifies
the county, and is still the county `NAME` for any sub-area. 

* **identity** the attribute value uniquely identifies the geometry as a whole, there are no other geometries with the same value

Arbitrary sub-areas will lose the **identity** property but becomes
a **constant** attribute. An example is:

* any point inside a county is still part of the county and must have the same value for county name, but it does not longer represent the (entire) geometry corresponding to that county.

We can specify the AGR of an attribute in an `sf` object by `st_set_agr`:
```{r}
nc <- system.file("gpkg/nc.gpkg", package="sf") %>%
	read_sf() %>%
	st_transform(32119)
nc1 <- nc %>% select(BIR74, SID74, NAME) %>%
	st_set_agr(c(BIR74 = "aggregate", SID74 = "aggregate", NAME = "identity"))
```
This helps to get rid of warnings that a particular attribute is assumed to be constant over a geometry, if it already is. The following no longer generates a warning
```{r}
nc1 %>% select(NAME) %>%
	st_centroid() %>%
	head(1)
```
and also changes AGR for `NAME` from `identity` to `constant` when replacing the geometry with the geometry's centroid:
```{r}
nc1 %>% select(NAME) %>%
	st_centroid() %>% 
	st_agr()
```

Identifying attribute-geometry relationships, and warning against
their absence is a first and simple implementation of the notion
that the types of phenomena we encounter in spatial data science
(like objects, fields, and aggregations) are not identified by their
geometrical representations (points, lines, polygons, rasters).
Making the wrong assumptions here easily leads to meaningless
analysis results [@stasch2014,@scheider2016].

## Spatial join

Spatial joins are similar to regular (left or inner) joins, where the
join criterion is not equality of one or more fields, but a spatial
predicate, such as that two records have intersecting geometries.
As an example, we can create a join between two tables,
```{r}
a = st_sf(a = 1:2, geom = st_sfc(st_point(c(0,0)), st_point(c(1,1))))
b = st_sf(b = 3:4, geom = st_sfc(st_linestring(rbind(c(2,0), c(0,2))), st_point(c(1,1))))
st_join(a, b)
st_join(a, b, left = FALSE)
st_join(b, a)
```
We see that unless `left = FALSE`, we get all elements (and
geometries) from the first argument, augmented with fields of
the second argument when geometries match. The example shows the
case where there are two geometries matching to point (1,1). The
spatial join predicate function can be freely chosen, e.g. from
the binary predicates listed in section \@ref(de9im).

When we match two sets of polygons, it may be a bit of a mess
to go through all the many matches. One way out of this is to only
provide the match with the largest overlap with the target geometry,
obtained by adding argument `largest = TRUE`. An example of this
is shown (visually) in figure \@ref(fig:largest).

(ref:foo) example of `st_join` with `largest = TRUE`: the label of the polygon in the top figure with the largest intersection with polygons in the bottom figure is assigned to the polygons of the bottom figure

```{r largest, out.width='60%', fig.cap='(ref:foo)', echo=FALSE}
# example of largest = TRUE:
nc <- st_transform(read_sf(system.file("shape/nc.shp", package="sf")), 2264)
gr = st_sf(
         label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "),
         geom = st_make_grid(nc))
gr$col = sf.colors(10, categorical = TRUE, alpha = .3)
# cut, to check, NA's work out:
gr = gr[-(1:30),]
suppressWarnings(nc_j <- st_join(nc, gr, largest = TRUE))
# the two datasets:
opar = par(mfrow = c(2,1), mar = rep(0,4))
plot(st_geometry(nc_j))
plot(st_geometry(gr), add = TRUE, col = gr$col)
text(st_coordinates(st_centroid(st_geometry(gr))), labels = gr$label)
# the joined dataset:
plot(st_geometry(nc_j), border = 'black', col = nc_j$col)
text(st_coordinates(st_centroid(st_geometry(nc_j))), labels = nc_j$label, cex = .8)
plot(st_geometry(gr), border = 'green', add = TRUE)
par(opar)
```

## Aggregate and summarise

Package `sf` provides `sf` methods for `stats::aggregate` and
`dplyr::summarise`. Both do essentially the same:

* given a grouping predicate (for `summarise`, obtained from `group_by`)
* given an aggregation function
* aggregate the selected attributes using this function, per group
* aggregate in addition the geometries.
* if `do_union` is `TRUE` (the default), union the aggregated geometries.

Unioning aggregated geometries dissolves for instance internal polygon
boundaries, which otherwise would lead to invalid `MULTIPOLYGON`
errors in subsequent analysis, or plotting of potentially unwanted
internal polygon boundaries. Figure \@ref(fig:union) illustrates
this.

```{r union, out.width='60%', fig.cap="left: invalid `MULTIPOLYGON` with two external rings with common boundary, right: valid `POLYGON` obtained after unioning the geometry left", echo=FALSE}
nc1 = st_geometry(nc)[1]
nc2 = st_geometry(nc)[2]
par(mfrow = c(1, 2))
plot(st_combine(c(nc1, nc2)))
plot(st_union(c(nc1, nc2)))
```

## Intersections

Suppose we have two datasets with different geometries and
attributes (left figure \@ref(fig:intersection)), and we want 
to compute their intersections:
```{r}
p1 = st_polygon(list(rbind(c(0,0), c(4,0), c(4,4), c(0,4), c(0,0))))
d1 = st_sf(a = c(3,1), geom = st_sfc(p1, p1 + c(4, 0)))
d2 = st_sf(b = c(4), geom = st_sfc(p1 * .75 + c(3, 2)))
```
What will the intersection of these two objects give?
```{r}
(i = st_intersection(d1, d2))
```
```{r intersection, out.width='60%', fig.cap="left: overlapping geometries (d2: red); right: intersection areas (i: grey)"}
plot(d1, xlim = c(0,8), ylim = c(0, 6), col = NA, border = 1, reset = FALSE)
plot(d2, col = NA, border = 'red', add = TRUE, lwd = 2)
plot(d1, xlim = c(0,8), ylim = c(0, 6), col = NA, border = 1, lwd = 2, reset = FALSE)
plot(d2, col = NA, border = 'red', add = TRUE, lwd = 3)
plot(st_geometry(i), add = TRUE, col = grey(c(.7,.9)), , border = 'green', lwd = 1)
```
As we see, this gives the areas of intersection, along with the
corresponding attributes for both contributing objects, and a
warning that attributes were assumed to be spatially constant.
Although this may be convenient in some cases, it may be entirely
meaningless in others. For instance, in case attribute `b` in object
`d2` represents the number of people living in `d2`, then after
the intersection we end up with twice as many people living over
a smaller area.

As seen in section \@ref(invalid), computing intersections easily
leads to errors caused by invalid geometries. Setting precision
(section \@ref(precision)) may prevent this.

## Area-weighted interpolation

Suppose we want to combine geometries and attributes of two datasets such that we get attribute values of the first dataset summarised for the geometries of the second. There are various ways we can go for this. The simplest one, building on the previous example, would be to obtain for the geometry of `d2` the attribute of `d1` that has the largest overlap with `d2`. This is obtained by
```{r}
st_join(d2, d1, largest = TRUE)
```

Another option would be to summarise the attribute, e.g. taking its mean, regardless the amount of overlap. This is obtained by
```{r}
aggregate(d1, d2, mean)
```

A third option is to apply area-weighted interpolation, meaning that we interpolate (average) the variable by taking into account the respective area contributions of overlap [@goodchild,@thomas,@do]. This is done e.g. by
```{r}
d3 = st_sfc(p1 * .75 + c(3, 2), p1 * .75 + c(3,3))
st_interpolate_aw(d1, d3, extensive = FALSE)$a
st_interpolate_aw(d1, d3, extensive = TRUE)$a
```

### Spatially intensive and extensive variables

The difference between the two examples for area-weighted interpolation is how the final weighted sum (value times area of intersection) is normalised: by the target area (extensive), or by the sum of the area covered (intensive, `extensive = FALSE`). Spatially intensive variables are variables for which the value, when we split an area, does not _principally_ change. An example might be temperature, elevation, or population density. Spatially extensive variables are variables for which the value is also split, according to the area. Examples are population (amount), or area.

## Exercises

* Add a variable to the `nc` dataset by `nc$State = "North Carolina"`. Which value should you attach to this variable for the attribute-geometry relationship (agr)?
* Create a new `sf` object from the geometry obtained by `st_union(nc)`, and assign `"North Carolina"` to the variable `State`. Which `agr` can you now assign to this attribute variable?
* Use `st_area` to add a variable with name `area` to `nc`. Compare the `area` and `AREA` variables in the `nc` dataset. What are the units of `AREA`? Are the two linearly related? If there are discrepancies, what could be the cause?
* Is the `area` variable intensive or extensive? Is its agr equal to `constant`, `identity` or `aggregate`?
* Find the name of the county that contains `POINT(-78.34046 35.017)`
* Find the names of all counties with boundaries that touch county `Sampson`.
* List the names of all counties that are less than 50 km away from county `Sampson`.