diff --git a/foundations.qmd b/foundations.qmd index e46a886..821d659 100644 --- a/foundations.qmd +++ b/foundations.qmd @@ -240,6 +240,7 @@ _ Manipulation: A demonstration using `dataRetrieval` package stream gage data from USGS: ```{r} +library(dataRetrieval) flows <- dataRetrieval::readNWISdv( siteNumbers = '14187200', parameterCd = "00060" @@ -721,6 +722,8 @@ summary(us_states['total_pop_15']) The geometry is **sticky**, which means that it is carried along by default! +You can read `sf` objects from shapefiles using `st_read()` or `read_sf()`. + ### Raster Data Model Raster data can be *continuous* (e.g. elevation, precipitation, atmospheric deposition) or it can be *categorical* (e.g. land use, soil type, geology type). Raster data can also be image based rasters which can be single-band or multi-band. You can also have a temporal component with space-time cubes. @@ -754,7 +757,7 @@ Try providing a different bounding box for an area of your choosing or changing ```{r raster-layer, message=FALSE, warning=FALSE, error=FALSE} library(terra) -myrast <- rast(ncol=10, nrow = 10, xmax=-116,xmin=-126,ymin=42,ymax=46) +myrast <- rast(ncol = 10, nrow = 10, xmax = -116, xmin = -126, ymin = 42, ymax = 46) myrast ``` ::: @@ -764,7 +767,7 @@ myrast ```{r} #| eval: false # just an example: -myrast <- rast(ncol=40, nrow = 40, xmax=-120,xmin=-136,ymin=52,ymax=56) +myrast <- rast(ncol=40, nrow = 40, xmax=-120, xmin=-136, ymin=52, ymax=56) ``` ::: @@ -976,7 +979,7 @@ A coordinate reference system for spatial data can be defined in different ways - 'authority:code' identifier + `EPSG:4326` -Proj4 used to be the standard crs identifier in R but the preferrable way is ty use the AUTHORITY:CODE method which `sf` as well as software like [QGIS](https://docs.qgis.org/3.16/en/docs/user_manual/working_with_projections/working_with_projections.html) will recognize. The AUTHORITY:CODE method is durable and easily discoverable online. +Proj4 used to be the standard crs identifier in R but the preferable way is to use the AUTHORITY:CODE method which `sf` as well as software like [QGIS](https://docs.qgis.org/3.16/en/docs/user_manual/working_with_projections/working_with_projections.html) will recognize. The AUTHORITY:CODE method is durable and easily discoverable online. The `WKT` representation of EPSG:4326 in `sf` is: ```{r} @@ -1004,8 +1007,8 @@ st_crs(gages_sf)==st_crs(pnw) # transform one to the other gages_sf <- st_transform(gages_sf, st_crs(pnw)) ggplot() + - geom_sf(data=gages_sf, color="blue") + - geom_sf(data=pnw, color="black", fill=NA) + + geom_sf(data = gages_sf, color="blue") + + geom_sf(data = pnw, color="black", fill=NA) + labs(title="USGS Stream Gages in the Pacific Northwest") + theme_bw() ``` @@ -1099,7 +1102,7 @@ See the [Geocomputation with R s2 section for further details](https://r.geocomp :::{.callout-note} There are sometimes good reasons for turning `S2` off during an R session. See issue 1771 in sf’s GitHub repo - the default behavior can make code that would work with S2 turned off (and with older versions of sf) fail. -Situations that can cause problems include operations on polygons that are not valid according to S2’s stricter definition. If you see error message such as #> Error in s2_geography_from_wkb ... you may want to turn `s2` off and try your operation again, perhaps turning of `s2` and running a topology clean operation like `st_make_valid()`. +Situations that can cause problems include operations on polygons that are not valid according to S2’s stricter definition. If you see error message such as `#> Error in s2_geography_from_wkb ...` you may want to turn `s2` off and try your operation again, perhaps turning of `s2` and running a topology clean operation like `st_make_valid()`. ::: ## Geographic Data I/O @@ -1157,11 +1160,11 @@ plot(citylims$geometry, axes=T, main='Oregon City Limits') # plot it! #### Solution `read_sf` is an `sf` alternative to st_read (see [this](https://keen-swartz-3146c4.netlify.com/intro.html#read) section 1.2.2). Try reading in `citylims` data above using `read_sf` and notice difference, and check out `help(read_sf)`. `read_sf` and `write_sf` are simply aliases for st_read and st_write with modified default arguments. Big differences are: -- stringsAsFactors=FALSE -- quiet=TRUE -- as_tibble=TRUE +- `stringsAsFactors=FALSE` +- `quiet=TRUE` +- `as_tibble=TRUE` -Note that stringsAsFactors = FALSE is the new default in R versions >= 4.0 +Note that `stringsAsFactors = FALSE` is the new default in R versions >= 4.0 #### Geodatabases We use `st_read` or `read_sf` similarly for reading in an ESRI file geodatabase feature: @@ -1285,8 +1288,8 @@ Load stock Landsat 7 .tif file that comes with package ```{r stars_read, message=FALSE, warning=FALSE, error=FALSE} library(stars) tif = system.file("tif/L7_ETMs.tif", package = "stars") -read_stars(tif) %>% - dplyr::slice(index = 1, along = "band") %>% +read_stars(tif) |> + dplyr::slice(index = 1, along = "band") |> plot() ``` @@ -1307,7 +1310,7 @@ In the steps below, we library(readr) library(ggplot2) library(Rspatialworkshop) -gages = read_csv(system.file("extdata/Gages_flowdata.csv", package = "Rspatialworkshop")) +gages = read.csv(system.file("extdata/Gages_flowdata.csv", package = "Rspatialworkshop")) gages_sf <- gages |> st_as_sf(coords = c("LON_SITE", "LAT_SITE"), crs = 4269, remove = FALSE) |> diff --git a/slides.html b/slides.html index a50b16c..89027d5 100644 --- a/slides.html +++ b/slides.html @@ -1,12 +1,12 @@ +!function(t,e){"object"==typeof exports&&"object"==typeof module?module.exports=e():"function"==typeof define&&define.amd?define([],e):"object"==typeof exports?exports.ClipboardJS=e():t.ClipboardJS=e()}(this,function(){return n={686:function(t,e,n){"use strict";n.d(e,{default:function(){return o}});var e=n(279),i=n.n(e),e=n(370),u=n.n(e),e=n(817),c=n.n(e);function a(t){try{return document.execCommand(t)}catch(t){return}}var f=function(t){t=c()(t);return a("cut"),t};var l=function(t){var e,n,o,r=1 - - + - + @@ -487,9 +486,9 @@ ul.task-list{list-style: none;} ul.task-list li input[type="checkbox"] { width: 0.8em; -margin: 0 0.8em 0.2em -1em; vertical-align: middle; +margin: 0 0.8em 0.2em -1.6em; +vertical-align: middle; } - pre > code.sourceCode { white-space: pre; position: relative; } pre > code.sourceCode > span { display: inline-block; line-height: 1.25; } pre > code.sourceCode > span:empty { height: 1.2em; } @@ -555,7 +554,7 @@ code span.vs { color: #20794d; } code span.wa { color: #5e5e5e; font-style: italic; } - + @@ -1710,7 +1708,7 @@ -

2023-10-18

+

10/18/23

Welcome!

  1. Please visit INSERT LINK HERE to view the workshop’s accompanying workbook

    2. @@ -1794,7 +1792,7 @@
    -
    +
    @@ -1822,7 +1820,8 @@

    -
    Figure 1: Distribution of Oregon counties.
    +

    Figure 1: Distribution of Oregon counties.

    +

Your Turn

@@ -1834,6 +1833,9 @@ mapview(my_town, col="red", col.regions = "red") + mapview(counties, alpha.regions = .1) + +
+

Spatial Data Structures

-

Figure 3: The simple features data model.

The sf package

+

Figure 3: The simple features data model.

The sf Package

@@ -1941,7 +1944,8 @@

-
Figure 5: sf LINESTRING
+

Figure 5: sf LINESTRING

+

@@ -1971,7 +1975,8 @@

-
Figure 6: sf POLYGON
+

Figure 6: sf POLYGON

+

Putting Them All Together

@@ -1986,7 +1991,8 @@

-
Figure 7: sf POINT, LINESTRING, and POLYGON overlain
+

Figure 7: sf POINT, LINESTRING, and POLYGON overlain

+

Spatial Data Frames

@@ -1994,6 +2000,8 @@
  • sf objects combine data.frame attributes with POINT, LINESTRING, and POLYGON building blocks (these are the geometry)
  • Also characterize dimensionality, bounding box, coordinate reference system, attribute headers
    • +
  • Read from shapefiles using st_read() or read_sf() +
    • library(spData)
@@ -2038,10 +2046,414 @@
 
    
 
    So far we have explored plotting using ggplot(), but you can also use sf’s base plotting via plot(). Read more about plotting in sf by running ?plot.sf and then make an appropriate plot of the us_states data.
    
 
    
+
    
+
    
+
    
+
    
+
+
    
+05:00
+
    
+
    
+
    
+

    Raster Data Model

    +
      +
    • Raster data can be continuous or categorical
      • +
    • Can be image based; have a temporal component
      • +
    + +

    Figure 9: The raster data model.

    The terra Package

    +
      +
    • We use the terra package to work with spatial raster data in R +
      • +
    • +terra is the modern replacement to raster +
      • +
    • Learn more about terra at their website https://rspatial.org/ +
      • +
    • +stars for spatio-temporal rasters
      • +

    Creating a Raster Object

    +
      +
    • Create an empty SpatRaster object; define rows, columns, bounding box
      • +
    +
    +
    library(terra)
+myrast <- rast(ncol = 10, nrow = 10, xmax = -116, xmin = -126, ymin = 42, ymax = 46)
+str(myrast)
    +
    +
    S4 class 'SpatRaster' [package "terra"]
    +
    +
    print(myrast)
    +
    +
    class       : SpatRaster 
+dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
+resolution  : 1, 0.4  (x, y)
+extent      : -126, -116, 42, 46  (xmin, xmax, ymin, ymax)
+coord. ref. : lon/lat WGS 84
    +
    +
    +

    Your Turn

    +
    +

    We just printed myrast and saw several components: class, dimensions, resolution, extent, and coord. ref. Explicitly return and inspect each of these pieces using class(), ncell(), res(), ext(), and crs(). Bonus: Return the number of raster cells using ncell().

    +
    +
    +
    +
    +
    + +
    +05:00 +
    +
    +
    +

    Manipulating Raster Objects

    +
      +
    • Let’s add values to the raster object
      • +
    +
    +
    values(myrast) <- 1:ncell(myrast)
    +
    +
    +
    plot(myrast)
    +
    + +

    Manipulating Raster Objects

    +
    +
    +

    +

    Figure 10: Raster plot.

    +
    +
    +
    +

    Reading a Raster Object

    +
      +
    • Lets read in a raster object from the spDataLarge package that contains elevation data from Zion National Park
      • +
    +
    +
    raster_filepath <- system.file("raster/srtm.tif", package = "spDataLarge")
+my_rast <- rast(raster_filepath)
    +
    +
    +
    plot(my_rast)
    +
    + +

    Reading a Raster Object

    +
    +
    +

    +

    Figure 11: Raster plot of elevation from Zion National Park.

    +
    +
    +
    +

    Your Turn

    +
    +

    The spatRaster object in terra can hold multiple layers. Read in the four bands of Landsat 8 data from Zion National Park by running

    +
    +
    landsat = system.file("raster/landsat.tif", package = "spDataLarge")
+landsat = rast(landsat)
    +
    +

    The plot this raster using plot().

    +
    +
    +
    +
    +
    + +
    +03:00 +
    +
    +
    +

    Spatio-Temporal Raster Objects

    +
      +
    • The stars package has tools for spatio-temporal raster data
      • +
    • +stars integrates with sf (recall that terra does not)
      • +
    • Learn more at their website https://r-spatial.github.io/stars/ +
      • +

    Coordinate Reference Systems

    +

    A coordinate reference system (CRS) is made up of several components

    +
      +
    1. Coordinate system: The x,y grid that defines where your data lies in space
      2. +
    2. Horizontal and vertical units: The units that describe the coordinate system
      3. +
    3. Datum: The modeled version of the earth’s shape
      4. +
    4. Projection details (if applicable): The mathematical equation used to flatten objects from round surface (earth) to flat surface (paper or screen)
      5. +
    +

    Properly specifying coordinate reference systems is a crucial step in handling or analyzing spatial data

    +

    The Ellipsoid and Geoid

    +
      +
    • The earth can be thought of as an ellipsoid defined by a semi-major (equatorial) axis and a semi-minor (polar) axis
      • +
    • More precisely, it is a geoid that is not perfectly smooth
      • +
    • A datum is built on top of the ellipsoid and incorporates local features
      • +
    • The most appropriate datum depends on the location of interest (in the US we typically use NAD83)
      • +

    The Ellipsoid and Geoid

    + +

    Figure 12: The ellipsoid and geoid.

    Your Turn

    +
    +

    There are several websites that are helpful for learning more about specific coordinate reference systems

    +
      +
    1. Spatial Reference: https://spatialreference.org/ +
      2. +
    2. Projection Wizard: https://projectionwizard.org/ +
      3. +
    3. EPSG: https://epsg.io/ +
      4. +
    +

    Spend a few minutes learning about one or more of these references by visiting their websites.

    +
    +
    +
    +
    +
    + +
    +04:00 +
    +
    +
    +

    Specifying a CRS

    +

    There are many ways to specify a CRS:

    +
      +
    1. Proj4 +
        +
      • +proj = longlat + ellps=WGS84 ...
        • +
      +
      2. +
    2. OGC WKT / ESRI WKT +
        +
      • GEOGCS["WGS 84", DATUM["WGS_1984, SPHERIOD["WGS 84 ...]]]
        • +
      +
      3. +
    3. authority:code identifier +
        +
      • EPSG: 4326
        • +
      +
      4. +
    +
      +
    • The authority:code identifier is the modern identifier in R, and our preferred method
      • +

    Projected Coordinate Systems

    +
      +
    • Projected coordinates have been projected to two-dimensional space according to a CRS
      • +
    • Projected coordinates have an origin, x-axis, y-axis, and unit of measure
      • +
    • Conformal projections preserve shape
      • +
    • Equal-area projections preserve area
      • +
    • Equidistant projections preserve distance
      • +
    • Azimuthal projection preserves direction
      • +

    Projected Coordinate Systems

    +
      +
    • Here an example using vector data; see the workbook for an example using raster data
      • +
    +
    +
    library(Rspatialworkshop)
+data("pnw")
+# transform one to the other
+utmz11 <- 2153
+pnw <- st_transform(pnw, crs = utmz11)
+ggplot() + 
+  geom_sf(data = pnw,  color="black", fill=NA) +
+  labs(title="State Boundaries in the Pacific Northwest") +
+  theme_bw()
    +
    + +

    Projected Coordinate Systems

    +
    +
    +

    +

    Figure 13: Gages and PNW data in the same projection.

    +
    +
    +
    +

    Your Turn

    +
    +

    Re-project the pnw data to a different projected CRS. Then plot using base R or ggplot2.

    +
    +
    +
    +
    +
    + +
    +06:00 +
    +
    +
    +

    A Note on S2 +

    +
      +
    • +sf version 1.0.0 supports spherical geometry operations via its interface to Google’s S2 spherical geometry engine
      • +
    • +S2 is an example of a Discrete Global Grid System (DGGS)
      • +
    • +sf can run with s2 on or off and by default the S2 geometry engine is turned on:
      • +
    +
    +
    sf::sf_use_s2()
    +
    +
    [1] TRUE
    +
    +
    +

    Geographic Data I/O (Input/Ouptut)

    +

    There are several ways to read spatial data into R

    +
      +
    1. Load spatial data from our machine or a remote source
      2. +
    2. Load spatial data as part of an R package
      3. +
    3. Load data using an API (which often makes use of an R package)
      4. +
    4. Convert flat files with x, y coordinates to spatial data
      5. +
    5. Geocoding data “by-hand” (we saw this earlier)
      6. +

    Vector Data I/O

    +

    sf can read numerous file types:

    +
      +
    1. shapefiles
      2. +
    2. geodatabases
      3. +
    3. geopackages
      4. +
    4. geojson
      5. +
    5. spatial database files
      6. +

    Read in a Shapefile

    +
    +
    filepath <- system.file("extdata/city_limits.shp", package = "Rspatialworkshop")
+citylims <- read_sf(filepath)
+plot(st_geometry(citylims), axes = T, main = "Oregon City Limits")
    +
    + +

    Read in a Shapefile

    +
    +
    +

    +

    Figure 14: Oregon City Limits

    +
    +
    +
    +

    Your Turn

    +
    +

    Run ?read_sf() and compare read_sf() to st_read(). Our preference is to use read_sf() – why do you think that is?

    +
    +
    +
    +
    +
    + +
    +03:00 +
    +
    +
    +

    Read in a Geodatabase

    +
    +
    filepath <- system.file("extdata/StateParkBoundaries.gdb", package = "Rspatialworkshop")
+# List all feature classes in a file geodatabase
+st_layers(filepath)
+# Read the feature class
+parks <- st_read(dsn = filepath, layer="StateParkBoundaries")
+ggplot(parks) + 
+  geom_sf()
    +
    + +

    Read in a Geodatabase

    +
    +
    +

    +

    Figure 15: State Park Boundaries

    +
    +
    +
    +

    Read in a Geopackage

    +
    +
    filepath <- system.file("gpkg/nc.gpkg", package="sf")
+nc <- read_sf(filepath)
+ggplot(nc) + 
+  geom_sf()
    +
    + +

    Read in a Geopackage

    +
    +
    +

    +

    Figure 16: North Carolina Counties

    +
    +
    +
    +

    Geopackage Advantages

    +

    Why use geopackages over shapefiles?

    +
      +
    1. Geopackages avoid mult-file format of shapefiles
      2. +
    2. Geopackages avoid the 2gb limit of shapefiles
      3. +
    3. Geopackages are open-source and follow OGC standards
      4. +
    4. Geopackages are lighter in file size than shapefiles
      5. +
    5. Geopackages avoid the 10-character limit to column headers in shapefile attribute tables (stored in archaic .dbf files)
      6. +

    Other Approaches

    +

    The workbook shows how to read in data using

    +
      +
    1. Open spatial data sources
      2. +
    2. +R packages
      3. +
    3. OpenStreetMap
      4. +

    Raster Data I/O

    +
      +
    • Read in data using rast() +
      • +
    +
    +
    filepath <- system.file("ex/elev.tif", package="terra")
+elev <- rast(filepath)
+barplot(elev, digits = -1, las = 2, ylab = "Frequency")
    +
    + +

    Raster Data I/O

    +
    +
    +

    +

    Figure 17: Elevation Data Barplot

    +
    +
    +
    +

    Raster Data I/O

    +
    +
    plot(elev)
    + +
    +

    Figure 18: Elevation Data

    Raster Data I/O

    +
    +
    library(stars)
+filepath <- system.file("tif/L7_ETMs.tif", package = "stars")
+L7 <- read_stars(filepath) |>
+  slice(index = 1, along = "band")
+plot(L7)
    +
    + +

    Raster Data I/O

    +
    +
    +

    +

    Figure 19: Landsat 7

    +
    +
    +
    +

    Convert Flat Files to Spatial

    +

    It is common for a flat file (e.g., a .csv) to have columns that indicate coordinates – how do we turn this into a spatial vector object?

    +
    +
    filepath <- system.file("extdata/Gages_flowdata.csv", package = "Rspatialworkshop")
+gages <- read.csv(filepath)
+gages_sf <- gages |>
+  st_as_sf(coords = c("LON_SITE", "LAT_SITE"), crs = 4269, remove = FALSE)
+ggplot(data = gages_sf) + 
+  geom_sf()
    +
    +

    Spatial Data Workshop

    -
    +

    Convert Flat Files to Spatial

    +
    +
    +

    +

    Figure 20: Stream Gages From a Flat File

    +
    +
    +
    +
    @@ -2508,9 +2920,6 @@ Reveal = _Reveal; install(); }; - Plugin.togglePdfExport = function () { - togglePdfExport(); - }; } return Plugin; @@ -2556,9 +2965,6 @@ downloadDrawings: revealMenuToolHandler(function () { RevealChalkboard.download(); }), - togglePdfExport: revealMenuToolHandler(function () { - PdfExport.togglePdfExport(); - }), };