class: center, middle, inverse, title-slide # Working with geospatial data ## Spatial computations ### Ben Baumer ### SDS 192April 13, 2020(
http://beanumber.github.io/sds192/lectures/mdsr_geo_07-computations.html
) --- class: center, middle, inverse # Spatial computations --- ## Area of polygons ```r # need lwgeom package! # install.packages("lwgeom") library(sf) ``` .footnote[https://github.com/r-spatial/lwgeom] -- ```r library(macleish) macleish_layers %>% pluck("buildings") %>% * st_area() ``` ``` ## Units: [m^2] ## [1] 38.35302 21.58409 102.00932 120.18765 87.40698 33.14448 44.25576 ## [8] 33.64470 163.47860 86.17566 115.16984 22.74444 207.11708 15.17230 ## [15] 55.09754 13.92451 24.31579 11.04709 248.12974 27.36442 292.06704 ## [22] 318.26350 195.68233 490.45395 315.18383 266.58585 307.19604 ``` --- ## Length of lines ```r macleish_layers %>% pluck("trails") %>% * st_length() ``` ``` ## Units: [m] ## [1] 699.63680 969.28430 899.85448 360.51939 831.63070 193.83152 ## [7] 187.50997 208.10487 1108.46789 66.91166 139.69878 69.13609 ## [13] 1347.04310 1228.17373 173.43296 ``` --- ## Distance between two objects ```r macleish_layers %>% pluck("challenge_courses") %>% * st_distance() ``` ``` ## Units: [m] ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] ## [1,] 0.00000 87.01874 58.16644 54.30083 40.61683 71.94224 79.03738 ## [2,] 87.01874 0.00000 64.40206 86.19723 53.37472 80.71333 112.32982 ## [3,] 58.16644 64.40206 0.00000 22.90191 25.59878 105.49286 125.97253 ## [4,] 54.30083 86.19723 22.90191 0.00000 38.65573 115.13225 130.50235 ## [5,] 40.61683 53.37472 25.59878 38.65573 0.00000 79.96865 100.78762 ## [6,] 71.94224 80.71333 105.49286 115.13225 79.96865 0.00000 33.21282 ## [7,] 79.03738 112.32982 125.97253 130.50235 100.78762 33.21282 0.00000 ## [8,] 134.63075 200.30901 192.68496 186.83130 171.97140 126.49255 93.86570 ## [9,] 96.41066 174.81594 153.72376 144.07113 136.62081 114.69075 88.62370 ## [10,] 91.56304 164.06823 149.72162 142.79317 130.31543 98.88768 70.96091 ## [11,] 67.70425 145.39951 125.60801 117.73497 107.48927 89.39660 68.16075 ## [12,] 56.91141 142.07067 111.44346 99.94277 97.36270 99.88291 85.62784 ## [,8] [,9] [,10] [,11] [,12] ## [1,] 134.63075 96.41066 91.56304 67.70425 56.91141 ## [2,] 200.30901 174.81594 164.06823 145.39951 142.07067 ## [3,] 192.68496 153.72376 149.72162 125.60801 111.44346 ## [4,] 186.83130 144.07113 142.79317 117.73497 99.94277 ## [5,] 171.97140 136.62081 130.31543 107.48927 97.36270 ## [6,] 126.49255 114.69075 98.88768 89.39660 99.88291 ## [7,] 93.86570 88.62370 70.96091 68.16075 85.62784 ## [8,] 0.00000 50.54635 44.34930 69.69217 92.37203 ## [9,] 50.54635 0.00000 19.01418 29.56191 44.82205 ## [10,] 44.34930 19.01418 0.00000 25.34287 48.91735 ## [11,] 69.69217 29.56191 25.34287 0.00000 25.50183 ## [12,] 92.37203 44.82205 48.91735 25.50183 0.00000 ``` --- ## Put a buffer around something ```r macleish_layers %>% pluck("streams") %>% * st_buffer(dist = 0.1) ``` ``` ## Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle = ## endCapStyle, : st_buffer does not correctly buffer longitude/latitude data ``` ``` ## dist is assumed to be in decimal degrees (arc_degrees). ``` ``` ## Simple feature collection with 13 features and 1 field ## geometry type: POLYGON ## dimension: XY ## bbox: xmin: -72.78641 ymin: 42.33478 xmax: -72.57133 ymax: 42.55891 ## CRS: EPSG:4326 ## First 10 features: ## Id geometry ## 1 1 POLYGON ((-72.74708 42.5258... ## 2 1 POLYGON ((-72.77373 42.4610... ## 3 1 POLYGON ((-72.76165 42.4026... ## 4 1 POLYGON ((-72.77642 42.4679... ## 5 1 POLYGON ((-72.765 42.40847,... ## 6 1 POLYGON ((-72.77245 42.4273... ## 7 1 POLYGON ((-72.77425 42.4584... ## 8 3 POLYGON ((-72.77909 42.4637... ## 9 3 POLYGON ((-72.71389 42.5449... ## 10 3 POLYGON ((-72.77938 42.4785... ``` --- ## Units!! - calculations are in the units of the projection! ```r macleish_layers %>% pluck("streams") %>% st_crs() ``` ``` ## Coordinate Reference System: ## User input: EPSG:4326 ## wkt: ## GEOGCS["WGS 84", ## DATUM["WGS_1984", ## SPHEROID["WGS 84",6378137,298.257223563, ## AUTHORITY["EPSG","7030"]], ## AUTHORITY["EPSG","6326"]], ## PRIMEM["Greenwich",0, ## AUTHORITY["EPSG","8901"]], ## UNIT["degree",0.0174532925199433, ## AUTHORITY["EPSG","9122"]], ## AUTHORITY["EPSG","4326"]] ``` - Important: Need interpretable units --- ## Workaround -
Step 1: Convert to Albers Equal Area - area-preserving - in meters - centered on US .footnote[https://epsg.io/102008] -- -
Convert back to EPSG:4326 for plotting ```r proj4_aea <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs" streams_wbuffer <- macleish_layers %>% pluck("streams") %>% st_transform(proj4_aea) %>% # put a 10 meter buffer around the streams st_buffer(dist = 10) %>% st_transform(4326) ``` --- ## Can use alternate tiles ```r library(leaflet) leaflet() %>% * addProviderTiles(providers$Esri.NatGeoWorldMap) %>% addPolylines(data = streams_wbuffer) ```
.footnote[http://leaflet-extras.github.io/leaflet-providers/preview/index.html] --- background-image: url(https://raw.githubusercontent.com/rstudio/cheatsheets/master/pngs/sf.png) background-size: contain .footnote[https://github.com/rstudio/cheatsheets/blob/master/sf.pdf]