Skip to content

Commit

Permalink
fixing typos, code formating and header fomating in spatial lesson
Browse files Browse the repository at this point in the history
  • Loading branch information
camilavargasp committed Mar 1, 2024
1 parent 0b35aee commit ac372b5
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 17 deletions.
46 changes: 29 additions & 17 deletions materials/sections/geospatial-vector-analysis.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -109,10 +109,10 @@ class(ak_regions)

`sf` objects usually have two types of classes: `sf` and `data.frame`.

Unlike a typical `data.frame`, an `sf` object has spatial metadata (`geometry type`, `dimension`, `bbox`, `epsg (SRID)`, `proj4string`) and an additional column typically named `geometry` that contains the spatial data.

Since our shapefile object has the `data.frame` class, viewing the contents of the object using the `head()` function or other exploratory functions shows similar results as if we read in data using `read.csv()` or `read_csv()`.

But, unlike a typical `data.frame`, an `sf` object has spatial metadata (`geometry type`, `dimension`, `bbox`, `epsg (SRID)`, `proj4string`) and an additional column typically named `geometry` that contains the spatial data.

```{r}
#| message: false
head(ak_regions)
Expand Down Expand Up @@ -242,7 +242,9 @@ In this case, we want to find what region each city falls within, so we will use

```{r}
#| eval: false
pop_joined <- st_join(pop_4326, ak_regions_3338, join = st_within)
pop_joined <- st_join(pop_4326,
ak_regions_3338,
join = st_within)
```

This gives an error!
Expand All @@ -254,11 +256,14 @@ Error: st_crs(x) == st_crs(y) is not TRUE
Turns out, this won't work right now because our coordinate reference systems are not the same. Luckily, this is easily resolved using `st_transform()`, and projecting our population object into Alaska Albers.

```{r}
pop_3338 <- st_transform(pop_4326, crs = 3338)
pop_3338 <- st_transform(pop_4326,
crs = 3338)
```

```{r}
pop_joined <- st_join(pop_3338, ak_regions_3338, join = st_within)
pop_joined <- st_join(pop_3338,
ak_regions_3338,
join = st_within)
head(pop_joined)
```
Expand All @@ -271,7 +276,7 @@ There are many different types of joins you can do with geospatial data. Examine

**3. Calculate the total population by region using `group_by()` and `summarize()`**

Next we compute the total population for each region. In this case, we want to do a `group_by()` and `summarise()` as this were a regular `data.frame`. Otherwise all of our point geometries would be included in the aggregation, which is not what we want. Our goal is just to get the total population by region. We remove the sticky geometry using `as.data.frame()`, on the advice of the `sf::tidyverse` help page.
Next we compute the total population for each region. In this case, we want to do a `group_by()` and `summarize()` as this were a regular `data.frame`. Otherwise all of our point geometries would be included in the aggregation, which is not what we want. Our goal is just to get the total population by region. We remove the sticky geometry using `as.data.frame()`, on the advice of the `sf::tidyverse` help page.

```{r}
pop_region <- pop_joined %>%
Expand All @@ -285,7 +290,9 @@ head(pop_region)
And use a regular `left_join()` to get the information back to the Alaska region shapefile. Note that we need this step in order to regain our region geometries so that we can make some maps.

```{r}
pop_region_3338 <- left_join(ak_regions_3338, pop_region, by = "region")
pop_region_3338 <- left_join(ak_regions_3338,
pop_region,
by = "region")
# plot to check
plot(pop_region_3338["total_pop"])
Expand Down Expand Up @@ -330,7 +337,7 @@ Save the spatial object to disk using `write_sf()` and specifying the filename.
write_sf(pop_region_3338, "data/ak_regions_population.shp")
```

### Visualize with `ggplot`
## Visualize with `ggplot`

`ggplot2` now has integrated functionality to plot sf objects using `geom_sf()`.

Expand All @@ -347,7 +354,10 @@ ggplot(pop_region_3338) +
theme_bw()
```

We can also plot multiple shapefiles in the same plot. Say if we want to visualize rivers in Alaska, in addition to the location of communities, since many communities in Alaska are on rivers. We can read in a rivers shapefile, doublecheck the `crs` to make sure it is what we need, and then plot all three shapefiles - the regional population (polygons), the locations of cities (points), and the rivers (linestrings).
We can also plot multiple shapefiles in the same plot. Say if we want to visualize rivers in Alaska, in addition to the location of communities, since many communities in Alaska are on rivers. We can read in a rivers shapefile, doublecheck the `crs` to make sure it is what we need, and then plot all three shapefiles
- the regional population (polygons),
- the locations of cities (points), and
- the rivers (linestrings).

```{r}
#| echo: false
Expand All @@ -359,19 +369,22 @@ st_crs(rivers_3338)

```{r}
#| eval: false
rivers_3338 <- read_sf(here("/data/shapefiles/ak_rivers_simp.shp"))
rivers_3338 <- read_sf("/data/shapefiles/ak_rivers_simp.shp")
st_crs(rivers_3338)
```

Note that although no EPSG code is set explicitly, with some sluething we can determine that this is `EPSG:3338`. [This site](https://epsg.io) is helpful for looking up EPSG codes.

```{r}
ggplot() +
geom_sf(data = pop_region_3338, aes(fill = total_pop)) +
geom_sf(data = pop_3338, size = 0.5) +
geom_sf(data = pop_region_3338,
aes(fill = total_pop)) +
geom_sf(data = pop_3338,
size = 0.5) +
geom_sf(data = rivers_3338,
aes(linewidth = StrOrder)) +
scale_linewidth(range = c(0.05, 0.5), guide = "none") +
scale_linewidth(range = c(0.05, 0.5),
guide = "none") +
labs(title = "Total Population by Alaska Region",
fill = "Total Population") +
scale_fill_continuous(low = "khaki",
Expand All @@ -381,7 +394,7 @@ ggplot() +
```


## Incorporate base maps into static maps using `ggspatial`
### Incorporate base maps into static maps using `ggspatial`

The `ggspatial` package has a function that can add tile layers from a few predefined tile sources like OpenStreetMap. The tiles will get projected into the CRS of the `sf` object pass into geom_sf(). Therefore no transformation is needed in this case.

Expand Down Expand Up @@ -423,8 +436,6 @@ ggplot(data = pop_3857) + # pop polygon layer
```




## Visualize `sf` objects with `leaflet`

We can also make an interactive map from our data above using `leaflet`.
Expand All @@ -451,7 +462,8 @@ st_crs(pop_region_3338)
Since `leaflet` requires that we use an unprojected coordinate system, let's use `st_transform()` yet again to get back to WGS84.

```{r}
pop_region_4326 <- pop_region_3338 %>% st_transform(crs = 4326)
pop_region_4326 <- pop_region_3338 %>%
st_transform(crs = 4326)
```

```{r}
Expand Down
2 changes: 2 additions & 0 deletions materials/session_19.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@ title-block-banner: true





{{< include /sections/geospatial-vector-analysis.qmd >}}


0 comments on commit ac372b5

Please sign in to comment.