From ac372b55a2fef30c0f1ed5dd508ea2cd592e3483 Mon Sep 17 00:00:00 2001 From: camilavargasp Date: Fri, 1 Mar 2024 08:04:30 -0800 Subject: [PATCH] fixing typos, code formating and header fomating in spatial lesson --- .../sections/geospatial-vector-analysis.qmd | 46 ++++++++++++------- materials/session_19.qmd | 2 + 2 files changed, 31 insertions(+), 17 deletions(-) diff --git a/materials/sections/geospatial-vector-analysis.qmd b/materials/sections/geospatial-vector-analysis.qmd index 602e2094..cae3583f 100644 --- a/materials/sections/geospatial-vector-analysis.qmd +++ b/materials/sections/geospatial-vector-analysis.qmd @@ -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) @@ -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! @@ -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) ``` @@ -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 %>% @@ -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"]) @@ -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()`. @@ -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 @@ -359,7 +369,7 @@ 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) ``` @@ -367,11 +377,14 @@ Note that although no EPSG code is set explicitly, with some sluething we can de ```{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", @@ -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. @@ -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`. @@ -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} diff --git a/materials/session_19.qmd b/materials/session_19.qmd index 25adff1d..70ac59f3 100644 --- a/materials/session_19.qmd +++ b/materials/session_19.qmd @@ -6,6 +6,8 @@ title-block-banner: true + + {{< include /sections/geospatial-vector-analysis.qmd >}}