A04-Vector-Operations.Rmd

# Vector operations in R{-}

```{r, child="_session-info.Rmd", echo = FALSE}
```

We'll first load spatial objects used in this exercise from a remote website: a polygon object that delineates  Maine counties; another polygon object that delineates distances to Augusta (Maine) as concentric circles; and a line object that shows the highway system that runs through Maine. These data are already stored as R data objects thus eliminating the need for any data conversion.

```{r}
z <- gzcon(url("http://colby.edu/~mgimond/Spatial/Data/Income_schooling.rds"))
s1 <- readRDS(z)

z <- gzcon(url("http://colby.edu/~mgimond/Spatial/Data/dist.rds"))
s2 <- readRDS(z)

z <- gzcon(url("http://colby.edu/~mgimond/Spatial/Data/highway.rds"))
l1 <- readRDS(z)
```

A map of the above layers is shown below.

```{r, out.width=450}
library(tmap)

# Get min/max bounding box for both polygon shapes
library(sp)
b1 <- bbox(s1)
b2 <- bbox(s2)
b3 <- pmax(b1,b2)
b3[,1] <- pmin(b1[,1],b2[,1])

# Plot all three features (note the bbox=b3 parameter that ensures that  
# the map extent encompasses all polygon extents)
tm_shape(s1, bbox = b3) + tm_fill(col="grey") + tm_borders(col = "white") +
  tm_shape(s2) + tm_fill(alpha = 0.2, col = "red") + tm_borders(col = "red") +
  tm_shape(l1) + tm_lines(col="yellow")
```

The attributes table for both polygon objects are shown next. Note that each shape object has a unique set of attributes as well as a unique number of records

```{r eval=T, echo=FALSE, results='asis', out.width=550, fig.cap="Attribute tables for the Maine spatial object, `s1`, (left table) and the distance to Augusta spatial object, `s2` (right table)."}
library(gridExtra)
tb.me  <- tableGrob( data.frame(s1),  theme = ttheme_default(base_size = 11))
tb.cir <- tableGrob( data.frame(s2),  theme = ttheme_default(base_size = 11))
tb <- gtable_combine(tb.me,tb.cir, along=1)
grid.arrange(tb)
```

## Dissolve by contiguous shape {-#app4_1}

To dissolve all polygons that share at least one line segment, simply pass the object name to `raster`'s `aggregate` function. In this example, we dissolve all polygons to create a single outline of the state of Maine.

```{r, out.width=400}
library(raster)
ME <- aggregate(s1)
qtm(ME)
```

 
## Dissolve by attribute {-#app4_2}

First, we'll create a new column whose value is binary (TRUE/FALSE) depending on whether or not the county income is below the counties' median income value.

```{r, out.width=400}
s1$med <- s1$Income > median(s1$Income)
tm_shape(s1) + tm_fill(col="med") +tm_borders(col = "white")
```

Next, we'll dissolve all polygons by the `med` value. Any polygons sharing at least one line segment having the same `med` value will be dissolved into a single polygon.

```{r, out.width=450}
ME.inc <- aggregate(s1, by= "med")
tm_shape(ME.inc) + tm_fill(col="med") +tm_borders()
```

The aggregation function will, by default, elliminate all other attribute values. If you wish to summarize other attribute values along with the attribute used for aggregation use `dplyr`'s piping operation. For example, to compute the median `Income` value for each of the below/above median income groups type the following:

```{r echo=FALSE}
# Unloading raster seems to be necessary for raster::intersect() to work
# properly later on in the script when running bookdown (the aggregate and intersect
# functions get masked by the base functions.)
detach("package:raster", unload=TRUE)
```

To summarize education by below/above median income

```{r, out.width=450}
library(dplyr)
ME.inc$Income <- s1@data %>% group_by(med) %>%
                          summarize(medinc = median(Income)) %>% .$medinc
tm_shape(ME.inc) + tm_fill(col="Income") +tm_borders()
```

To view the attributes table with both the aggregate variable, `med`, and the median income variable, `Income`, type:

```{r}
ME.inc@data
```


## Calculate area {-#app4_3}

To calculate a polygon's area, you can use `rgeos`'s `gArea` function. For example, to compute the area in km^2^, type the following:

```{r}
library(rgeos)
ME.inc$Area <- gArea(ME.inc, byid=TRUE) / 1000000
```

The `gArea` function computes the area in map units which happens to be in meters in our example. To convert the area from m^2^ to km^2^, we simply divide by `1,000,000` in the above chunk of code. This produces the following attributes table.

```{r echo=FALSE}
# Remove unloading rgeos seems to be necessary for raster::intersect() to work
# properly when running bookdown
detach("package:rgeos", unload=TRUE)
```

 
```{r}
ME.inc@data
```

## Subsetting {-#app4_4}

You can use conventional R dataframe manipulation operations to subset by attribute values. For example, to subset by county name (e.g. `Kennebec` county), type:

```{r, out.width=400}
ME.ken <- s1[s1$NAME == "Kennebec",]
qtm(ME.ken)
```

To subset by a range of attribute values (e.g. subset by income values that are less than the median value), type:

```{r, out.width=400}
ME.inc2 <- s1[s1$Income < median(s1$Income), ]
qtm(ME.inc2)
```

## Intersecting {-#app4_5}

To intersect two polygon objects, you can use `raster`'s `intersect` function.

```{r, out.width=400}
library(raster)
clp1 <- intersect(s1,s2)
tm_shape(clp1) + tm_fill(col="Income") +tm_borders()
```

`intersect` keeps all features that overlap along with their combined attributes. Note that new polygons are created which can increase the size of the attributes table beyond the size of the combined input attributes table.

```{r tblgrb2, echo=FALSE, results='asis', out.width=500, fig.height=9.5}
library(gridExtra)
tb1 <- tableGrob(data.frame(clp1), rows = NULL, theme=ttheme_default(base_size = 10))
grid.arrange(tb1)
```

### Intersecting polygons with line and point objects {-}

This method will also intersect line and point objects with a polygon object, but the output will be either a line/point or a polygon depending on the order in which the objects are passed to `intersect`. If the line object is passed first, the output will be a line object that falls within the polygons of the input polygon object. If the polygon object is passed first, then all polygons that have a line segment within their boundaries will be returned. For example:

To output all counties having at least one segment of the line feature running through them, type:

```{r out.width=450}
library(raster)
clp2 <- intersect(s1,l1)
tm_shape(clp2) + tm_fill(col="grey") +tm_borders() + 
  tm_shape(l1) + tm_lines(col="red")
```

To output all line segments that fall within the concentric distance circles of `s2`, type:

```{r out.width=450}
library(raster)
clp3 <- intersect(l1,s2)
tm_shape(s2) + tm_fill(col="grey") +tm_borders() + 
  tm_shape(clp3) + tm_lines(col="red")
```

In both cases, the output shape objects inherit both the original attribute values as well as the attributes from the intersecting object.

## Unioning {-#app4_6}

To union two polygon objects, use `raster`'s `union` function. For example,

```{r, out.width=450}
library(raster)
un1  <- union(s1,s2)
tm_shape(un1) + tm_fill(col="Income") + tm_borders()
```

This produces the following attributes table.
 
```{r echo=FALSE, results='asis', out.width=500, fig.height=13}
library(gridExtra)
tb2 <- tableGrob(data.frame(un1), rows = NULL, theme=ttheme_default(base_size = 9))
grid.arrange(tb2)
```