class: center, middle, inverse, title-slide # Spatial data visualization ## Statistical Computing & Programming ### Shawn Santo --- ## Supplementary materials Full video lecture available in Zoom Cloud Recordings Additional resources - Simple Features for R [vignettes](https://r-spatial.github.io/sf/) - [CRS in R](https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf) by Melanie Frazier - [Leaflet for R](https://rstudio.github.io/leaflet/) --- ## Data Data for today is available at - Package `sf` - NC counties, births, sudden infant deaths ``` nc <- st_read(system.file("shape/nc.shp", package = "sf")) ``` - [NC OneMap](https://www.nconemap.gov) - lots of spatial data on all things NC - [NC PMTW Streams 2020](https://www.nconemap.gov/datasets/0cc135e6c6244c9e9646b45ee3cb6c1e_0) - [NC game lands](https://www.nconemap.gov/datasets/ncwrc::game-lands-general) - [London cholera data](http://blog.rtwilson.com/john-snows-cholera-data-in-more-formats/) --- class: inverse, center, middle # A little history --- ## Spatial data analysis then John Snow's map that changed the world... <br/> <center> <img src="images/cholera_map.jpg"> </center> Backstory: read more [here](https://www.theguardian.com/news/datablog/2013/mar/15/john-snow-cholera-map) --- ## Spatial data analysis now Spatial data analysis and visualization was a critical component in Operation Neptune Spear. <center> <img src="images/bin_laden_compound.jpeg" width=600 height=400> </center> Backstory: read more [here](https://www.theatlantic.com/politics/archive/2011/05/the-little-known-agency-that-helped-kill-bin-laden/238454/) --- class: inverse, center, middle # Introduction --- ## Spatial data is different Our **typical tidy data frame**: .tiny[ ``` #> # A tibble: 336,776 x 19 #> year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time #> <int> <int> <int> <int> <int> <dbl> <int> <int> #> 1 2013 1 1 517 515 2 830 819 #> 2 2013 1 1 533 529 4 850 830 #> 3 2013 1 1 542 540 2 923 850 #> 4 2013 1 1 544 545 -1 1004 1022 #> 5 2013 1 1 554 600 -6 812 837 #> 6 2013 1 1 554 558 -4 740 728 #> 7 2013 1 1 555 600 -5 913 854 #> 8 2013 1 1 557 600 -3 709 723 #> 9 2013 1 1 557 600 -3 838 846 #> 10 2013 1 1 558 600 -2 753 745 #> # … with 336,766 more rows, and 11 more variables: arr_delay <dbl>, #> # carrier <chr>, flight <int>, tailnum <chr>, origin <chr>, dest <chr>, #> # air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm> ``` ] --- ## Spatial data is different A **simple features object**: .tiny[ ``` #> Simple feature collection with 100 features and 5 fields #> geometry type: MULTIPOLYGON #> dimension: XY #> bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965 #> geographic CRS: NAD27 #> First 10 features: #> AREA PERIMETER CNTY_ CNTY_ID NAME geometry #> 1 0.114 1.442 1825 1825 Ashe MULTIPOLYGON (((-81.47276 3... #> 2 0.061 1.231 1827 1827 Alleghany MULTIPOLYGON (((-81.23989 3... #> 3 0.143 1.630 1828 1828 Surry MULTIPOLYGON (((-80.45634 3... #> 4 0.070 2.968 1831 1831 Currituck MULTIPOLYGON (((-76.00897 3... #> 5 0.153 2.206 1832 1832 Northampton MULTIPOLYGON (((-77.21767 3... #> 6 0.097 1.670 1833 1833 Hertford MULTIPOLYGON (((-76.74506 3... #> 7 0.062 1.547 1834 1834 Camden MULTIPOLYGON (((-76.00897 3... #> 8 0.091 1.284 1835 1835 Gates MULTIPOLYGON (((-76.56251 3... #> 9 0.118 1.421 1836 1836 Warren MULTIPOLYGON (((-78.30876 3... #> 10 0.124 1.428 1837 1837 Stokes MULTIPOLYGON (((-80.02567 3... ``` ] --- Another **simple features** object: .tiny[ ``` #> Simple feature collection with 1858 features and 2 fields #> geometry type: MULTILINESTRING #> dimension: XYZ #> bbox: xmin: -9383722 ymin: 4162088 xmax: -8947168 ymax: 4381832 #> z_range: zmin: 0 zmax: 0 #> projected CRS: WGS 84 / Pseudo-Mercator #> First 10 features: #> Displ_Name FIRST_WRC_ #> 1 Big Sandy Creek Delayed Harvest Trout Waters #> 2 Helton Creek Delayed Harvest Trout Waters #> 3 Little River Delayed Harvest Trout Waters #> 4 Cranberry Creek Hatchery Supported Trout Waters #> 5 Big Pine Creek Hatchery Supported Trout Waters #> 6 Bledsoe Creek Hatchery Supported Trout Waters #> 7 Brush Creek Hatchery Supported Trout Waters #> 8 <NA> Hatchery Supported Trout Waters #> 9 (Big) Glade Creek Hatchery Supported Trout Waters #> 10 Little River Hatchery Supported Trout Waters #> geometry #> 1 MULTILINESTRING Z ((-901773... #> 2 MULTILINESTRING Z ((-907349... #> 3 MULTILINESTRING Z ((-903290... #> 4 MULTILINESTRING Z ((-904529... #> 5 MULTILINESTRING Z ((-901483... #> 6 MULTILINESTRING Z ((-903607... #> 7 MULTILINESTRING Z ((-901813... #> 8 MULTILINESTRING Z ((-902513... #> 9 MULTILINESTRING Z ((-902642... #> 10 MULTILINESTRING Z ((-903124... ``` ] --- Yet another **simple features** object: .tiny[ ``` #> Simple feature collection with 94 features and 2 fields #> geometry type: MULTIPOLYGON #> dimension: XY #> bbox: xmin: 127456.7 ymin: 26544.91 xmax: 923528.7 ymax: 318097.4 #> projected CRS: NAD83 / North Carolina #> First 10 features: #> GML_HAB SUM_ACRES geometry #> 1 Alcoa 11109.5589 MULTIPOLYGON (((512096.2 18... #> 2 Alligator River 24439.0891 MULTIPOLYGON (((869633.1 24... #> 3 Angola Bay 34067.3821 MULTIPOLYGON (((713079.4 11... #> 4 Bachelor Bay 2786.2577 MULTIPOLYGON (((813742.2 23... #> 5 Bertie County 3881.4660 MULTIPOLYGON (((797133.8 24... #> 6 Bladen Lakes State Forest 33671.8439 MULTIPOLYGON (((658970.6 95... #> 7 Brinkleyville 3189.7974 MULTIPOLYGON (((714741 2769... #> 8 Buckhorn 491.3477 MULTIPOLYGON (((589723.7 25... #> 9 Buckridge 17965.7187 MULTIPOLYGON (((871137.4 21... #> 10 Buffalo Cove 6630.9453 MULTIPOLYGON (((381445.9 26... ``` ] --- ## Spatial data plotting needs care <img src="lec_11_files/figure-html/unnamed-chunk-6-1.png" width="100%" style="display: block; margin: auto;" /> ``` #> null device #> 1 ``` --- <img src="lec_11_files/figure-html/unnamed-chunk-8-1.png" width="100%" style="display: block; margin: auto;" /> --- <img src="lec_11_files/figure-html/unnamed-chunk-9-1.png" width="100%" style="display: block; margin: auto;" /> --- class: center, middle ## Can we combine the two plots? --- ## What happened? <img src="lec_11_files/figure-html/unnamed-chunk-10-1.png" width="100%" style="display: block; margin: auto;" /> --- class: center, middle ## We can combine them, but more care is needed --- <img src="lec_11_files/figure-html/unnamed-chunk-11-1.png" width="100%" style="display: block; margin: auto;" /> --- ## Spatial data challenges 1. Different data types exist. 2. Special attention must be given to the coordinate reference system (CRS). 3. Manipulating spatial data objects is similar but not identical to manipulating data frame objects. --- class: inverse, center, middle # Spatial data and R --- ## Analysis of spatial data in R .pull-left[ <br/> - Package `raster` contains classes and tools for handling spatial raster data. <br/><br/> - Package `sf` combines the functionality of `sp`, `rgdal`, and `rgeos` into a single package based on tidy simple features. ] .pull-right[ ![](images/vector_raster_comparison.png) ] <br/> Whether or not you use vector or raster data depends on the type of problem and the data source. Our focus will be on vector data and package `sf`. *Source:* https://commons.wikimedia.org/wiki/File:Raster_vector_tikz.png --- ## Installing package `sf` From https://r-spatial.github.io/sf/index.html **Windows** Installing `sf` from source works under windows when Rtools is installed. This downloads the system requirements from rwinlib. **MacOS** ```bash brew install pkg-config brew install gdal ``` Once gdal is installed, you will be able to install sf package from source in R. **Linux** For Unix-alikes, GDAL (>= 2.0.1), GEOS (>= 3.4.0) and Proj.4 (>= 4.8.0) are required. --- ## Features and simple features - A **feature** is a thing or object in the real world: a house, a city, a park, a forest, etc. <br/><br/> - A **simple feature** as defined by OpenGIS Abstract specification is to have both spatial and non-spatial attributes. Spatial attributes are geometry valued, and simple features are based on 2D geometry with linear interpolation between vertices. <br/><br/> .tiny[ ```r Simple feature collection with 1858 features and 2 fields geometry type: MULTILINESTRING dimension: XYZ bbox: xmin: -9383722 ymin: 4162088 xmax: -8947168 ymax: 4381832 z_range: zmin: 0 zmax: 0 projected CRS: WGS 84 / Pseudo-Mercator First 10 features: Displ_Name FIRST_WRC_ geometry *1 Big Sandy Creek Delayed Harvest Trout Waters MULTILINESTRING Z ((-901773... 2 Helton Creek Delayed Harvest Trout Waters MULTILINESTRING Z ((-907349... 3 Little River Delayed Harvest Trout Waters MULTILINESTRING Z ((-903290... 4 Cranberry Creek Hatchery Supported Trout Waters MULTILINESTRING Z ((-904529... 5 Big Pine Creek Hatchery Supported Trout Waters MULTILINESTRING Z ((-901483... 6 Bledsoe Creek Hatchery Supported Trout Waters MULTILINESTRING Z ((-903607... 7 Brush Creek Hatchery Supported Trout Waters MULTILINESTRING Z ((-901813... 8 <NA> Hatchery Supported Trout Waters MULTILINESTRING Z ((-902513... 9 (Big) Glade Creek Hatchery Supported Trout Waters MULTILINESTRING Z ((-902642... 10 Little River Hatchery Supported Trout Waters MULTILINESTRING Z ((-903124... ``` ] --- ## Simple features examples <img src="lec_11_files/figure-html/unnamed-chunk-14-1.png" width="100%" style="display: block; margin: auto;" /> ``` #> null device #> 1 ``` --- ## `sf` objects .tiny[ ```r nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE) nc ``` ``` #> Simple feature collection with 100 features and 14 fields #> geometry type: MULTIPOLYGON #> dimension: XY #> bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965 #> geographic CRS: NAD27 #> First 10 features: #> AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 #> 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1 #> 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0 #> 3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5 #> 4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1 #> 5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9 #> 6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7 #> 7 0.062 1.547 1834 1834 Camden 37029 37029 15 286 0 #> 8 0.091 1.284 1835 1835 Gates 37073 37073 37 420 0 #> 9 0.118 1.421 1836 1836 Warren 37185 37185 93 968 4 #> 10 0.124 1.428 1837 1837 Stokes 37169 37169 85 1612 1 #> NWBIR74 BIR79 SID79 NWBIR79 geometry #> 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3... #> 2 10 542 3 12 MULTIPOLYGON (((-81.23989 3... #> 3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3... #> 4 123 830 2 145 MULTIPOLYGON (((-76.00897 3... #> 5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3... #> 6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3... #> 7 115 350 2 139 MULTIPOLYGON (((-76.00897 3... #> 8 254 594 2 371 MULTIPOLYGON (((-76.56251 3... #> 9 748 1190 2 844 MULTIPOLYGON (((-78.30876 3... #> 10 160 2038 5 176 MULTIPOLYGON (((-80.02567 3... ``` ] --- ## Class and other attributes: `sf` ```r class(nc) ``` ``` #> [1] "sf" "data.frame" ``` ```r names(attributes(nc)) ``` ``` #> [1] "names" "row.names" "class" "sf_column" "agr" ``` --- ## `sfc` objects ```r nc_polygons <- st_geometry(nc) nc_polygons ``` ``` #> Geometry set for 100 features #> geometry type: MULTIPOLYGON #> dimension: XY #> bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965 #> geographic CRS: NAD27 #> First 5 geometries: ``` ``` #> MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3... ``` ``` #> MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3... ``` ``` #> MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3... ``` ``` #> MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36... ``` ``` #> MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3... ``` --- ## Class and other attributes: `sfc` ```r class(nc_polygons) ``` ``` #> [1] "sfc_MULTIPOLYGON" "sfc" ``` ```r names(attributes(nc_polygons)) ``` ``` #> [1] "n_empty" "crs" "class" "precision" "bbox" ``` <br/> -- We see that `nc` has a class attribute `sf`, and object `nc_polygons` has a class attribute `sfc`. What methods are available? --- ```r sloop::s3_methods_class("sf") %>% print(n = 20) ``` ``` #> # A tibble: 92 x 4 #> generic class visible source #> <chr> <chr> <lgl> <chr> #> 1 [ sf FALSE registered S3method #> 2 [[<- sf FALSE registered S3method #> 3 $<- sf FALSE registered S3method #> 4 aggregate sf FALSE registered S3method #> 5 anti_join sf FALSE registered S3method #> 6 arrange sf FALSE registered S3method #> 7 as.data.frame sf FALSE registered S3method #> 8 cbind sf FALSE registered S3method #> 9 distinct sf FALSE registered S3method #> 10 dplyr_reconstruct sf FALSE registered S3method #> 11 filter sf FALSE registered S3method #> 12 full_join sf FALSE registered S3method #> 13 gather sf FALSE registered S3method #> 14 group_by sf FALSE registered S3method #> 15 group_split sf FALSE registered S3method #> 16 identify sf FALSE registered S3method #> 17 inner_join sf FALSE registered S3method #> 18 left_join sf FALSE registered S3method #> 19 merge sf FALSE registered S3method #> 20 mutate sf FALSE registered S3method #> # … with 72 more rows ``` --- ```r sloop::s3_methods_class("sfc") %>% print(n = 20) ``` ``` #> # A tibble: 65 x 4 #> generic class visible source #> <chr> <chr> <lgl> <chr> #> 1 [ sfc FALSE registered S3method #> 2 [<- sfc FALSE registered S3method #> 3 as.data.frame sfc FALSE registered S3method #> 4 c sfc FALSE registered S3method #> 5 format sfc FALSE registered S3method #> 6 fortify sfc FALSE registered S3method #> 7 identify sfc FALSE registered S3method #> 8 obj_sum sfc FALSE registered S3method #> 9 Ops sfc FALSE registered S3method #> 10 print sfc FALSE registered S3method #> 11 rep sfc FALSE registered S3method #> 12 scale_type sfc FALSE registered S3method #> 13 st_area sfc FALSE registered S3method #> 14 st_as_binary sfc FALSE registered S3method #> 15 st_as_grob sfc FALSE registered S3method #> 16 st_as_s2 sfc FALSE registered S3method #> 17 st_as_sf sfc FALSE registered S3method #> 18 st_as_text sfc FALSE registered S3method #> 19 st_bbox sfc FALSE registered S3method #> 20 st_boundary sfc FALSE registered S3method #> # … with 45 more rows ``` --- ## Reading and writing spatial data - `st_read()` / `st_write()`, Shapefile, GeoJSON, KML, ... - `st_as_sfc()` - `st_as_text()`, well-known text format - `st_as_binary()`, well-known binary format <br/><br/><br> See https://r-spatial.github.io/sf/articles/sf2.html for the full set of driver availability. --- ## Plotting with `plot()` ```r plot(nc, max.plot = 14) ``` <img src="lec_11_files/figure-html/unnamed-chunk-24-1.png" style="display: block; margin: auto;" /> --- ```r plot(nc["NAME"]) ``` <img src="lec_11_files/figure-html/unnamed-chunk-25-1.png" style="display: block; margin: auto;" /> --- ```r par(oma=c(0,2,0,0)) plot(nc["AREA"], graticule = TRUE, axes = TRUE, las = 1) ``` <img src="lec_11_files/figure-html/unnamed-chunk-26-1.png" style="display: block; margin: auto;" /> --- ## What is happening with `[` and the `sf` object? ```r nc["AREA"] ``` ``` #> Simple feature collection with 100 features and 1 field #> geometry type: MULTIPOLYGON #> dimension: XY #> bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965 #> geographic CRS: NAD27 #> First 10 features: #> AREA geometry #> 1 0.114 MULTIPOLYGON (((-81.47276 3... #> 2 0.061 MULTIPOLYGON (((-81.23989 3... #> 3 0.143 MULTIPOLYGON (((-80.45634 3... #> 4 0.070 MULTIPOLYGON (((-76.00897 3... #> 5 0.153 MULTIPOLYGON (((-77.21767 3... #> 6 0.097 MULTIPOLYGON (((-76.74506 3... #> 7 0.062 MULTIPOLYGON (((-76.00897 3... #> 8 0.091 MULTIPOLYGON (((-76.56251 3... #> 9 0.118 MULTIPOLYGON (((-78.30876 3... #> 10 0.124 MULTIPOLYGON (((-80.02567 3... ``` --- ```r nc$AREA ``` ``` #> [1] 0.114 0.061 0.143 0.070 0.153 0.097 0.062 0.091 0.118 0.124 0.114 0.153 #> [13] 0.143 0.109 0.072 0.190 0.053 0.199 0.081 0.063 0.044 0.064 0.086 0.128 #> [25] 0.108 0.170 0.111 0.180 0.104 0.077 0.142 0.059 0.131 0.122 0.080 0.118 #> [37] 0.219 0.118 0.155 0.069 0.066 0.145 0.134 0.100 0.099 0.116 0.201 0.180 #> [49] 0.094 0.134 0.168 0.106 0.168 0.207 0.144 0.094 0.203 0.141 0.070 0.065 #> [61] 0.146 0.142 0.154 0.118 0.078 0.125 0.181 0.143 0.091 0.130 0.103 0.095 #> [73] 0.078 0.104 0.098 0.091 0.060 0.131 0.241 0.082 0.120 0.172 0.121 0.163 #> [85] 0.138 0.098 0.167 0.204 0.121 0.051 0.177 0.080 0.195 0.240 0.125 0.225 #> [97] 0.214 0.240 0.042 0.212 ``` --- ```r par(oma=c(0,2,0,0)) plot(nc["AREA"], col = "lightblue", graticule = TRUE, axes = TRUE, las = 1) ``` <img src="lec_11_files/figure-html/unnamed-chunk-29-1.png" style="display: block; margin: auto;" /> --- ```r par(oma=c(0,2,0,0)) plot(st_geometry(nc), graticule = TRUE, axes = TRUE, las = 1) ``` <img src="lec_11_files/figure-html/unnamed-chunk-30-1.png" style="display: block; margin: auto;" /> --- ## Plotting with `ggplot()` ```r ggplot(nc) + * geom_sf() + theme_bw(base_size = 16) ``` <img src="lec_11_files/figure-html/unnamed-chunk-31-1.png" style="display: block; margin: auto;" /> --- ```r ggplot(nc) + * geom_sf(aes(fill = AREA)) + scale_fill_gradient(low = "#fee8c8", high = "#7f0000") + theme_bw(base_size = 16) ``` <img src="lec_11_files/figure-html/unnamed-chunk-32-1.png" style="display: block; margin: auto;" /> --- ## Patchwork works <img src="lec_11_files/figure-html/unnamed-chunk-33-1.png" style="display: block; margin: auto;" /> --- ## Patchwork works Aside: What is wrong with the last set of figures? -- <br/> Code for previous slide: ```r p1 <- ggplot(nc) + geom_sf(aes(fill = SID74)) + scale_fill_gradient(low = "#fff7f3", high = "#49006a") + theme_bw(base_size = 16) p2 <- ggplot(nc) + geom_sf(aes(fill = SID79)) + scale_fill_gradient(low = "#fff7f3", high = "#49006a") + theme_bw(base_size = 16) p1 / p2 ``` --- ## Plotting with `leaflet` ```r library(leaflet) ``` -- ```r # create a color palette function pal <- colorNumeric(palette = "Blues", domain = nc$BIR74) leaflet(nc, width = "100%") %>% addTiles() %>% addPolygons(color = "grey") %>% addPolygons(stroke = FALSE, fillOpacity = 0.8, smoothFactor = 0.2, color = ~pal(BIR74)) %>% addLegend(position = "bottomright", pal = pal, values = ~BIR74, opacity = 0.8) ``` --- ## Plotting with `leaflet`
--- ## Exercise Use `ggplot` to create a choropleth map for the proportion of sudden infant deaths, for the period of July 1, 1974 to June 30, 1979. <img src="lec_11_files/figure-html/unnamed-chunk-38-1.png" width="100%" style="display: block; margin: auto;" /> ??? ```r nc %>% select(BIR74, SID74) %>% mutate(SID74_prop = SID74 / (BIR74 + SID74)) %>% ggplot() + geom_sf(aes(fill = SID74_prop)) + scale_fill_gradient(low = "#fff7f3", high = "#49006a") + labs(title = "July 1, 1974 to June 30, 1979", fill = "", subtitle = "Proportion of SID by county") + theme_void() + theme(plot.margin = margin(0, 1, 0, 1, "cm")) ``` --- class: inverse, center, middle # Map layers --- ## PMTW streams data The North Carolina Wildlife Resources Commission developed the Public Mountain Trout Waters (PMTW) digital data to enhance planning and management of trout waters. The GIS dataset depicts the trout regulations in effect on trout waters (streams and impoundments) managed under the PMTW program as listed in 2020-2021. We can directly download and unzip the shapefile via ```r download.file(str_c("https://opendata.arcgis.com/datasets/", "0cc135e6c6244c9e9646b45ee3cb6c1e_0.zip", "?outSR=%7B%22latestWkid%22%3A3857%2C%22", "wkid%22%3A102100%7D"), destfile = "data/pmtw.zip") unzip("data/pmtw.zip", exdir = "data/") ``` -- To see the available files ```r list.files(path = "data/", pattern = "PMTW_*") ``` ``` #> [1] "PMTW_Streams_2020.cpg" "PMTW_Streams_2020.dbf" "PMTW_Streams_2020.prj" #> [4] "PMTW_Streams_2020.shp" "PMTW_Streams_2020.shx" "PMTW_Streams_2020.xml" ``` --- ## Read in the shapefile ```r nc_pmtw <- st_read("data/PMTW_Streams_2020.shp", quiet = TRUE) ``` ```r glimpse(nc_pmtw) ``` ``` #> Rows: 1,858 #> Columns: 11 #> $ FID <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,… #> $ Excel_ID <int> 0, 0, 2, 4, 5, 6, 7, 8, 8, 9, 10, 11, 12, 12, 13, 20, 23, … #> $ Displ_Name <chr> "Big Sandy Creek", "Helton Creek", "Little River", "Cranbe… #> $ MHTW_Reach <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA… #> $ FIRST_Reg_ <chr> "Helton Creek", "Helton Creek", "Little River", "Cranberry… #> $ FIRST_Reg1 <chr> "S.R. 1372 bridge to North Fork New River", "S.R. 1372 bri… #> $ FIRST_WRC_ <chr> "Delayed Harvest Trout Waters", "Delayed Harvest Trout Wat… #> $ FIRST_MHTW <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA… #> $ WEB_Refere <chr> "http://www.ncwildlife.org/Fishing/LawsSafety/FishingRegul… #> $ Shape__Len <dbl> 37.44214, 17.63478, 4996.62215, 2925.42681, 5796.97361, 11… #> $ geometry <MULTILINESTRING [m]> MULTILINESTRING Z ((-901773..., MULTILINES… ``` --- ## Metadata for each `sf` object `nc`: ```r Simple feature collection with 100 features and 14 fields geometry type: MULTIPOLYGON dimension: XY bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965 *geographic CRS: NAD27 ``` `nc_pmtw`: ```r Simple feature collection with 1858 features and 10 fields geometry type: MULTILINESTRING dimension: XYZ bbox: xmin: -9383722 ymin: 4162088 xmax: -8947168 ymax: 4381832 z_range: zmin: 0 zmax: 0 *projected CRS: WGS 84 / Pseudo-Mercator ``` --- ## Check the CRS ```r st_crs(nc) ``` ```r Coordinate Reference System: User input: NAD27 wkt: GEOGCRS["NAD27", DATUM["North American Datum 1927", ELLIPSOID["Clarke 1866",6378206.4,294.978698213898, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["latitude",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["longitude",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4267]] ``` --- .tiny[ ```r st_crs(nc_pmtw) ``` ```r Coordinate Reference System: User input: WGS 84 / Pseudo-Mercator wkt: PROJCRS["WGS 84 / Pseudo-Mercator", BASEGEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]], CONVERSION["Popular Visualisation Pseudo-Mercator", METHOD["Popular Visualisation Pseudo Mercator", ID["EPSG",1024]], PARAMETER["Latitude of natural origin",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["False easting",0, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",0, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["easting (X)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing (Y)",north, ORDER[2], LENGTHUNIT["metre",1]], USAGE[ SCOPE["unknown"], AREA["World - 85°S to 85°N"], BBOX[-85.06,-180,85.06,180]], ID["EPSG",3857]] ``` ] --- ## Coordinate reference systems (CRS) - CRS provide a standardized way of describing locations. - Different CRS arise from various ways data were gathered, the locations, and purposes of the data. - A CRS is comprised of - an ellipsoid, to define the earth's shape; - a datum, to define the origin and orientation of coordinate axes; - a projection, to go from 3D to 2D. - It is important that you transform your spatial data to a common CRS before plotting. --- ## Transform CRS ```r nc_pmtw <- st_transform(nc_pmtw, crs = st_crs(nc)) ``` -- Check they are equal: ```r st_crs(nc) == st_crs(nc_pmtw) ``` ``` #> [1] TRUE ``` --- class: center, middle <img src="lec_11_files/figure-html/unnamed-chunk-52-1.png" style="display: block; margin: auto;" /> --- ## Code: map overlay with `plot()` ```r plot(st_geometry(nc), axes = T, las = 1, main = "NC PMTW Streams", cex.main = 3, cex.lab = 2, cex.axis = 1.5) plot(st_geometry(nc_pmtw), * add = T, col = "#3BB3D0") legend("bottomleft", legend = "Trout streams", fill = "#3BB3D0") ``` --- class: center, middle <img src="lec_11_files/figure-html/unnamed-chunk-54-1.png" style="display: block; margin: auto;" /> --- ## Code: map overlay with `ggplot2` ```r # use some fancy fonts library(extrafont) loadfonts() ggplot(nc) + geom_sf(alpha = 0.2) + geom_sf(data = nc_pmtw, aes(color = FIRST_WRC_), size = 1.5) + labs(title = "NC PMTW Streams", color = "Stream designation", caption = str_c("Public mountain trout waters are only in ", "western North Carolina")) + theme_void() + xlim(c(-84.5, -80.5)) + # zoom long ylim(c(34.8, 36.8)) + # zoom lat # custom theme theme(legend.position = c(0.25, 0.80), plot.title = element_text(size = 50, family = "Impact"), legend.text = element_text(size = 16, family = "Comic Sans MS"), legend.title = element_text(size = 20, family = "Comic Sans MS"), plot.caption = element_text(size = 12, family = "Comic Sans MS")) ``` --- ## Exercise Recreate John Snow's cholera map showing the cholera death locations and the water pump locations. Download the data (you should have a directory named `data/`): ```r download.file("http://rtwilson.com/downloads/SnowGIS_SHP.zip", destfile = "data/john_snow.zip") unzip("data/john_snow.zip", exdir = "data/") ``` Read in the data: ```r cholera <- st_read("data/SnowGIS_SHP/Cholera_Deaths.shp", quiet = TRUE) pumps <- st_read("data/SnowGIS_SHP/Pumps.shp", quiet = TRUE) ``` Before you start, check that both `sf` objects have the same CRS. --- ## References 1. Interactive Viewing of Spatial Data in R. (2021). https://r-spatial.github.io/mapview/index.html. 2. "Leaflet For R - Introduction". Rstudio.Github.Io, 2021, https://rstudio.github.io/leaflet/. 3. Melanie Frazier. Coordinate Reference Systems in R. https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf. 4. Simple Features for R. (2021). https://r-spatial.github.io/sf/.