class: center, middle, inverse, title-slide # Introduction to Geospatial Analysis in R ### Thomas Roh
Data Scientist, HDR ### 2018-09-05 --- class: inverse ## Software Requirements - R >= 3.4.3 - R Packages - sf >= 0.6-3 - sp >= 1.3-1 - ggplot2 >= 3.0.0.9 - leaflet >= 2.0.1 - raster >= 2.6-7 - httr >= 1.3.1 - rvest - RStudio Deskop (preferred) - Any other IDE for R will work as well - [Census API Key](https://api.census.gov/data/key_signup.html) --- class: inverse ## Outline - Introduction to R and GIS - GIS Data - Shapefiles, GeoJSON, WKT, rasters - CRS - spatial geometries - Reading Spatial Data into R - Data Representations - sp Package - sf Package - raster Package - leaflet Package - ggplot2 Package --- ## Data Types in R -- ```r typeof(1L) ``` ``` [1] "integer" ``` -- ```r typeof(1.1) ``` ``` [1] "double" ``` -- ```r typeof(TRUE) ``` ``` [1] "logical" ``` -- ```r typeof("string") ``` ``` [1] "character" ``` --- ## Data Types in R -- ```r typeof(factor(c("string1", "string2"))) ``` ``` [1] "integer" ``` -- ```r factor(c("string1", "string2")) ``` ``` [1] string1 string2 Levels: string1 string2 ``` -- ```r ordered(c("string1", "string2")) ``` ``` [1] string1 string2 Levels: string1 < string2 ``` --- ## Data Structures in R -- ```r ## all must be of the same type c(1L, 2L) ``` ``` [1] 1 2 ``` -- ```r matrix(rnorm(4), ncol = 2) ``` ``` [,1] [,2] [1,] -1.267763 -0.4673907 [2,] 1.734005 -0.2678249 ``` -- ```r str(array(rnorm(8), dim = c(2, 2, 2))) ``` ``` num [1:2, 1:2, 1:2] -1.378 0.351 -1.386 -0.329 -1.006 ... ``` -- ```r ## different data types can be stored together sapply(list(1L, "string", FALSE, 1.1), typeof) ``` ``` [1] "integer" "character" "logical" "double" ``` --- ## Dataframes in R -- ```r ## special class of list df <- data.frame(x = 1:2, y = c("string1", c("string2"))) df ``` ``` x y 1 1 string1 2 2 string2 ``` -- ```r df[1,"y", drop = FALSE] ``` ``` y 1 string1 ``` -- ```r df[df$x == 1, names(df) == 'x'] ``` ``` [1] 1 ``` --- ## Subsetting in R -- ```r df['y'] ``` ``` y 1 string1 2 string2 ``` -- ```r df$y ``` ``` [1] string1 string2 Levels: string1 string2 ``` -- ```r all.equal(df$y, df[['y']]) ``` ``` [1] TRUE ``` -- [Advanced R - Data Structures](http://adv-r.had.co.nz/Data-structures.html) --- ## GIS Data Files (Shapefile) -- - Vector geospatial data format specification developed by ESRI to minimize storage space, improve shape drawing speed, and for interoperatibily with GIS software -- - Typically in distributed in a zip file and stored in a directory with at least the following files: + Main file (.shp) + Index file (.shx) + dBASE table (.dbf) -- - The other most significant file is the projection format (.prj) which supplies the coordinate reference system -- - There are many other [files](https://en.wikipedia.org/wiki/Shapefile#Overview) that may be in the directory, which can be ignored for now. -- [Source: ESRI](http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf) --- ## GIS Data Files (GeoJSON) -- - Spatial data format developed in line with JavaScript Object Notation (JSON). GeoJSON objects are standardized JSON objects. -- - All information about the is contained in a single file (.geojson) - FeatureCollection->Feature->Geometry - Geometry Types: Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon, and GeometryCollection - CRS information (optional) - JSON is a named list of object data, can be nested and assymetrical -- ```javascript { "type": "Feature", "geometry": { "type": "Point", "coordinates": [ [-96, 43] ] } } ``` -- [Source: IETF](https://tools.ietf.org/html/rfc7946) --- ## GIS Data Files (WKT) -- - Well-known Text (WKT) is a fully open source data specification for spatial objects that is designed to be human and machine readable. -- - A simple text markup language that is easy to store in a text file or in a database as a character string. For databases that do not support spatial objects, you can store as WKT as character strings and use an interpreter in R to read it. -- - The `print` method for geometry objects in the [sf](https://cran.r-project.org/web/packages/sf/index.html) package will show WKT representation in the console. -- ``` POINT (-95.9345 41.2565) LINESTRING (-96.05645 41.22345, -95.90417 41.11618) POLYGON ((-95.9345 41.2565, -95.8345 41.2565, -95.9345 41.3565, -95.9345 41.2565)) ``` -- [Source: Open Geospatial Consortium](http://docs.opengeospatial.org/is/12-063r5/12-063r5.html) --- ## Spatial Data Types -- ```r library(sf) set.seed(82) p <- st_point(c(-95.9345, 41.2565)) mp <- st_multipoint(matrix(rnorm(4, sd = .1) + rep(p, each = 2), ncol = 2)) l <- st_linestring(mp) ml <- st_multilinestring(list(l + rnorm(4, sd = .05), l - rnorm(4, sd = .05))) plgn <- st_polygon(list(matrix(rep(p, 4) + c(0,0,.1,0,0,.1,0,0), ncol = 2, byrow = TRUE))) mplgn <- st_multipolygon(list(plgn -.15 ,plgn -.2)) ``` --- ## Spatial Data Types ```r par(mar = rep(0, 4)) plot(p, xlim = p[1]+c(-.1, .01), ylim = p[2]+c(-.2, .1)) plot(mp, add = TRUE, col = 'red') plot(l, add = TRUE, col = 'blue') plot(ml, add = TRUE, col = 'orange') plot(plgn, add = TRUE, col = 'purple') plot(mplgn, add = TRUE, col = 'yellow') ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-20-1.png" style="display: block; margin: auto;" /> --- ## Spatial Data Types ```r par(mar = rep(0, 4)) gcol <- st_geometrycollection(list(p, mp, l, ml, plgn, mplgn)) plot(gcol, xlim = p[1]+c(-.1, .01), ylim = p[2]+c(-.2, .1)) ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-21-1.png" style="display: block; margin: auto;" /> --- ## Reading a Shapefile in with `sf` -- ```r states <- st_read('dataWS/states/cb_2017_us_state_500k.shp', stringsAsFactors = FALSE, quiet = TRUE) %>% st_transform(crs = 4326) ``` -- ```r str(states) ``` --- ## Changing the Coordinate Reference System ```r states <- st_transform(states, crs = 4326) st_crs(states) ``` ``` Coordinate Reference System: EPSG: 4326 proj4string: "+proj=longlat +datum=WGS84 +no_defs" ``` -- ```r states %>% st_transform(crs = 4326) %>% st_crs() ``` -- - Common CRS + [4326 (WGS 84)](https://epsg.io/4326) + [3857 (Web Mercator)](https://epsg.io/3857) + [4269 (NAD83)](https://epsg.io/4269) --- ## Mapping with ggplot2 ```r library(ggplot2) ggplot(states[!states$STUSPS %in% c('HI', 'AK', 'AS', 'GU', 'MP', 'PR', 'VI'), ]) + geom_sf(fill = 'transparent', color = 'blue', size = 1) + theme_bw() ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-26-1.png" style="display: block; margin: auto;" /> --- ## Omaha City Limits ```r cityLimits <- st_read('dataWS/city_limits.csv', stringsAsFactors = FALSE, quiet = TRUE, crs = 4326) bb <- st_bbox(cityLimits, crs = st_crs(4326)) tracts <- st_read('dataWS/census_tracts_31/cb_2017_31_tract_500k.shp', stringsAsFactors = FALSE, quiet = TRUE) %>% st_transform(crs = 4326) ``` --- ## Omaha City Limits ```r ggplot(cityLimits) + geom_sf(aes(fill = TOWN)) + geom_sf(data = tracts, fill = 'transparent') + coord_sf(xlim = c(bb$xmin, bb$xmax), ylim = c(bb$ymin, bb$ymax)) + scale_fill_viridis_d(option = 'B') + theme_bw() ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-28-1.png" style="display: block; margin: auto;" /> --- ## Filtering Spatial Data -- ```r omaha <- cityLimits[cityLimits$TOWN == 'Omaha',] tracts <- tracts[st_intersects(tracts, omaha, sparse = FALSE), ] ggplot(omaha) + geom_sf(aes(fill = TOWN)) + geom_sf(data = tracts, fill = 'transparent') + coord_sf(xlim = c(bb$xmin, bb$xmax), ylim = c(bb$ymin, bb$ymax)) + scale_fill_viridis_d(option = 'B') + theme_bw() ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-29-1.png" style="display: block; margin: auto;" /> --- ## U.S. Census API -- ```r library(httr) query_acs <- function(var, stateNo, apiKey = Sys.getenv("census_token")) { endpoint <- "https://api.census.gov/data/2016/acs/acs5" varString <- paste0(c("NAME", var), collapse = ",") apiCall <- paste0(endpoint, "?get=", varString, "&for=tract:*&in=state:", stateNo, "&key=", apiKey) response <- GET(apiCall) census <- content(response) census <- lapply(census[-1], function(x) { setNames(data.frame(x, stringsAsFactors = FALSE), nm = census[[1]]) }) census <- Reduce(rbind, census) census[['geo_id']] <- paste0(census$state, census$county, census$tract) census } census <- query_acs(c("B00001_001E", "B06011_001E"), 31) census[1,] ``` ``` NAME B00001_001E B06011_001E state 1 Census Tract 9654, Adams County, Nebraska 924 32048 31 county tract geo_id 1 001 965400 31001965400 ``` <b>US Census: </b> [User Guide](https://www.census.gov/content/dam/Census/data/developers/api-user-guide/api-guide.pdf), [Census Variables](https://api.census.gov/data/2016/acs/acs5/variables.html), [Examples](https://api.census.gov/data/2013/acs1/examples.html) --- ## Chloropleths -- ```r names(census)[1:3] <- c('description', 'population', 'median_income') census[['population']] <- as.numeric(census[['population']]) census[['median_income']] <- as.numeric(census[['median_income']]) tracts <- merge(tracts, census, by.x = 'GEOID', by.y = 'geo_id') tracts[['density']] <- tracts[['population']]/(tracts[['ALAND']]/1000000) ``` --- ```r ggplot(tracts) + geom_sf(aes(fill = density)) + scale_fill_viridis_c(name = 'Population Density (per km2)') + theme_bw() ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-32-1.png" style="display: block; margin: auto;" /> ```r tracts[['description']][which.max(tracts$density)] ``` ``` [1] "Census Tract 39, Douglas County, Nebraska" ``` --- ## Interactive Mapping ```r library(leaflet) tracts <- tracts[tracts$median_income > 0,] centroids <- st_centroid(st_union(omaha)) tractBaseMap <- leaflet(tracts) %>% addTiles() %>% setView(centroids[[1]][1], centroids[[1]][2], zoom = 10) tractBaseMap ```
[leaflet](https://rstudio.github.io/leaflet/) --- ## Color Palettes ```r pal <- colorNumeric("viridis", tracts$median_income) palBin <- colorBin('viridis', tracts$median_income, 5) palQuant <- colorQuantile('viridis', tracts$median_income, 5) tractMap <- tractBaseMap %>% addPolygons(color = ~pal(median_income), group = 'Continuous', weight = 1, fillOpacity = .5) %>% addLegend(pal = pal, values = ~median_income, group = 'Continuous Legend', position = 'bottomright') ``` --- ```r tractMap ```
--- ## Legends ```r tractMap <- tractMap %>% addPolygons(color = ~palBin(median_income), group = 'Bin', weight = 1, fillOpacity = .5) %>% addPolygons(color = ~palQuant(median_income), group = 'Quantile', weight = 1, fillOpacity = .5) %>% addLegend(pal = palBin, values = ~median_income, group = 'Bin Legend', position = 'bottomright') %>% addLegend(pal = palQuant, values = ~median_income, group = 'Quantile Legend', position = 'bottomleft') %>% addLayersControl(baseGroups = c('Continuous', 'Bin', 'Quantile'), overlayGroups = c('Continuous Legend', 'Bin Legend', 'Quantile Legend')) ``` --- ```r tractMap ```
--- ## Omaha Building Violations ```r violations <- st_read( 'dataWS/OPBCV/Omaha_Planning_Building_Code_Violations.shp', quiet = TRUE, stringsAsFactors = FALSE) %>% st_transform(crs = 4326) vColors <- ifelse(grepl('Closed', violations$StatusCurr), 'black', 'red') tdown <- makeAwesomeIcon('thumbs-down', library = 'fa', markerColor = 'white', iconColor = vColors, squareMarker = TRUE) ``` [Font Awesome](https://fontawesome.com) --- ## Markers/Points ```r tractMap <- tractBaseMap %>% addPolygons(color = ~pal(median_income), group = 'Continuous', weight = 1, fillOpacity = .5) %>% addLegend(pal = pal, values = ~median_income, group = 'Continuous Legend', position = 'bottomright') %>% addAwesomeMarkers(data = violations, icon = tdown, clusterOptions = markerClusterOptions(), popup = ~paste( "<b>Permit:</b>",PermitNum, "<br>", "<b>Type:</b>", CaseType, "<br>")) ``` --- ```r tractMap ```
--- ## Spatial Aggregation of Violations ```r violWithin <- st_within(violations, tracts, sparse = FALSE) nViol <- colSums(violWithin) nViol[1:5] ``` ``` [1] 64 83 66 47 71 ``` ```r sum(nViol) ``` ``` [1] 4325 ``` ```r nrow(violations) ``` ``` [1] 4331 ``` --- ```r tracts <- cbind(tracts, nViol) ggplot(tracts) + geom_sf(aes(fill = nViol)) + geom_sf(data = violations[ rowSums(violWithin) == 0, ], color = 'red') + scale_fill_viridis_c(option = 'plasma') + theme_bw() ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-42-1.png" style="display: block; margin: auto;" /> --- ## Year Structures Were Built ```r yrBuilt <- query_acs(paste0("B25034_", sprintf("%03.f", 1:11), "E"), 31) head(yrBuilt, 5) ``` ``` NAME B25034_001E B25034_002E 1 Census Tract 9654, Adams County, Nebraska 1866 12 2 Census Tract 9655, Adams County, Nebraska 1689 0 3 Census Tract 9656, Adams County, Nebraska 2246 5 4 Census Tract 9657, Adams County, Nebraska 728 0 5 Census Tract 9658, Adams County, Nebraska 1036 0 B25034_003E B25034_004E B25034_005E B25034_006E B25034_007E B25034_008E 1 27 338 282 241 333 140 2 8 140 36 77 301 173 3 23 131 354 52 361 221 4 0 15 41 14 47 42 5 0 16 0 4 75 128 B25034_009E B25034_010E B25034_011E state county tract geo_id 1 73 58 362 31 001 965400 31001965400 2 240 143 571 31 001 965500 31001965500 3 174 366 559 31 001 965600 31001965600 4 54 46 469 31 001 965700 31001965700 5 332 121 360 31 001 965800 31001965800 ``` ```r yrBuilt <- yrBuilt[ ,!names(yrBuilt) %in% c('state', 'county', 'tract', 'NAME')] ``` --- ## Querying Census Data Labels ```r query_acs_var <- function(var) { endpoint <- "https://api.census.gov/data/2016/acs/acs5/variables/" x <- lapply(paste0(endpoint, var), function(x) {content(GET(x))}) setNames(sapply(x, getElement, 'label'), var) } varNames <- query_acs_var(paste0("B25034_", sprintf("%03.f", 1:11), "E")) varNames[1:5] ``` ``` B25034_001E "Estimate!!Total" B25034_002E "Estimate!!Total!!Built 2014 or later" B25034_003E "Estimate!!Total!!Built 2010 to 2013" B25034_004E "Estimate!!Total!!Built 2000 to 2009" B25034_005E "Estimate!!Total!!Built 1990 to 1999" ``` --- ## Prepare the Data ```r varNames <- gsub(" ", "_", gsub('[A-Za-z]+\\!\\!', "", varNames)) varNames ``` ``` B25034_001E B25034_002E B25034_003E "Total" "Built_2014_or_later" "Built_2010_to_2013" B25034_004E B25034_005E B25034_006E "Built_2000_to_2009" "Built_1990_to_1999" "Built_1980_to_1989" B25034_007E B25034_008E B25034_009E "Built_1970_to_1979" "Built_1960_to_1969" "Built_1950_to_1959" B25034_010E B25034_011E "Built_1940_to_1949" "Built_1939_or_earlier" ``` ```r names(yrBuilt)[1:11] <- varNames[names(yrBuilt)[1:11]] yrBuilt[, -12] <- lapply(yrBuilt[, -12], as.numeric) yrBuilt <- na.omit(yrBuilt) tracts <- st_as_sf(tracts, crs = 4326) tracts <- merge(tracts, yrBuilt, by.x = 'GEOID', by = 'geo_id') tracts <- as.data.frame(tracts) class(tracts) ``` ``` [1] "data.frame" ``` --- ```r trainDT <- tracts[ , c('Built_1939_or_earlier', 'Built_1940_to_1949', 'Built_1950_to_1959', 'Built_1960_to_1969', 'Built_1970_to_1979', 'Built_1980_to_1989', 'Built_1990_to_1999', 'Built_2000_to_2009', 'Built_2010_to_2013', 'Built_2014_or_later', 'Total', 'density', 'median_income', 'nViol')] trainDT[['built_1949']] <- trainDT$Built_1939_or_earlier + trainDT$Built_1940_to_1949 trainDT[['built_1969']] <- trainDT$Built_1950_to_1959 + trainDT$Built_1960_to_1969 trainDT[['built_1989']] <- trainDT$Built_1970_to_1979 + trainDT$Built_1980_to_1989 trainDT[['built_1990']] <- trainDT$Built_1990_to_1999 + trainDT$Built_2000_to_2009 + trainDT$Built_2010_to_2013 + trainDT$Built_2014_or_later trainDT <- trainDT[-c(1:10)] ``` --- ## Train a Model ```r poissonModel <- glm(nViol ~ ., family = 'poisson', data = trainDT) summary(poissonModel) ``` ``` Call: glm(formula = nViol ~ ., family = "poisson", data = trainDT) Deviance Residuals: Min 1Q Median 3Q Max -9.1355 -1.6920 0.2232 1.2533 6.4402 Coefficients: (1 not defined because of singularities) Estimate Std. Error z value Pr(>|z|) (Intercept) 4.526e+00 6.822e-02 66.339 < 2e-16 *** Total -1.830e-04 6.183e-05 -2.960 0.00308 ** density -1.188e-03 1.856e-04 -6.398 1.57e-10 *** median_income -6.522e-05 2.261e-06 -28.848 < 2e-16 *** built_1949 1.254e-03 6.716e-05 18.680 < 2e-16 *** built_1969 5.714e-04 6.993e-05 8.171 3.05e-16 *** built_1989 4.775e-04 8.706e-05 5.484 4.15e-08 *** built_1990 NA NA NA NA --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 4104.8 on 161 degrees of freedom Residual deviance: 1061.4 on 155 degrees of freedom AIC: 1776.6 Number of Fisher Scoring iterations: 5 ``` --- ## Singularity ```r alias(poissonModel) ``` ``` Model : nViol ~ Total + density + median_income + built_1949 + built_1969 + built_1989 + built_1990 Complete : (Intercept) Total density median_income built_1949 built_1969 built_1990 0 1 0 0 -1 -1 built_1989 built_1990 -1 ``` --- ## Results ```r poissonModel <- glm(nViol ~ ., family = 'poisson', data = trainDT[,-1]) summary(poissonModel) ``` ``` Call: glm(formula = nViol ~ ., family = "poisson", data = trainDT[, -1]) Deviance Residuals: Min 1Q Median 3Q Max -9.1355 -1.6920 0.2232 1.2533 6.4402 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 4.526e+00 6.822e-02 66.339 < 2e-16 *** density -1.188e-03 1.856e-04 -6.398 1.57e-10 *** median_income -6.522e-05 2.261e-06 -28.848 < 2e-16 *** built_1949 1.071e-03 3.782e-05 28.332 < 2e-16 *** built_1969 3.884e-04 5.340e-05 7.273 3.51e-13 *** built_1989 2.945e-04 5.688e-05 5.177 2.25e-07 *** built_1990 -1.830e-04 6.183e-05 -2.960 0.00308 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 4104.8 on 161 degrees of freedom Residual deviance: 1061.4 on 155 degrees of freedom AIC: 1776.6 Number of Fisher Scoring iterations: 5 ``` --- ## Model Diagnostics ```r par(mfrow = c(2,2)) plot(poissonModel) ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-50-1.png" style="display: block; margin: auto;" /> --- ## Model Error ```r tracts[['pred_nViol']] <- predict(poissonModel, type = "response") tracts[['pred_error']] <- tracts$nViol-tracts$pred_nViol qplot(poissonModel$y, predict(poissonModel, type = "response")) + geom_abline(slope = 1, intercept = 0, color = 'red') + labs(x = 'Actual', y = 'Predicted') + theme_bw() ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-51-1.png" style="display: block; margin: auto;" /> --- ## Spatial Comparison .pull-left[ ```r ggplot(tracts) + geom_sf(aes(fill = nViol)) + labs(fill = 'Violations') + scale_fill_viridis_c() + theme_bw() ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-52-1.png" style="display: block; margin: auto;" /> ] .pull-right[ ```r ggplot(tracts) + geom_sf(aes(fill = pred_nViol)) + labs(fill = 'Violations') + scale_fill_viridis_c() + theme_bw() ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-53-1.png" style="display: block; margin: auto;" /> ] --- ## Spatial Error ```r ggplot(tracts) + geom_sf(aes(fill = pred_error)) + scale_fill_viridis_c(option = 'cividis') + theme_bw() ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-54-1.png" style="display: block; margin: auto;" /> --- ## GIS Data Files (Rasters) -- - Rasters graphics are a very generic object class -- - For GIS analysis, a spatial raster is divided into equal units of the CRS -- - A raster is an array with pixel location on an X-Y plane and has topological information contained at the each location -- ```r aNames <- c(lapply(c("X", "Y"), paste0, 1:2), list(c("R", "G", "B"))) array(rnorm(10), dim = c(2,2,3), dimnames = aNames) ``` ``` , , R Y1 Y2 X1 -0.03789851 -0.9039875 X2 -0.31576897 0.4371677 , , G Y1 Y2 X1 -0.3617482 1.200864 X2 0.6089799 0.646111 , , B Y1 Y2 X1 -0.2890454 -0.03789851 X2 1.5248314 -0.31576897 ``` --- ## Working with Rasters -- - [raster](https://cran.r-project.org/web/packages/raster/raster.pdf) -- - sf does not have raster abilities yet -- - Works on large files because it can divide and conquer -- - RasterLayer, RasterStack, RasterBrick -- - [rasterfile](https://cran.r-project.org/web/packages/raster/vignettes/rasterfile.pdf) --- ## Bounding Box for Omaha Area ```r library(raster) coords <- centroids[[1]] bb <- c(xmin = -96.268204, xmax = -95.742829, ymin = 41.109533, ymax = 41.393063) plot(extent(bb)) ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-56-1.png" style="display: block; margin: auto;" /> --- ## Shuttle Radar Topography Mission ```r srtm <- getData('SRTM', lat = coords[2], lon = coords[1], path = 'dataWS/altitude', download = TRUE) srtm <- crop(srtm, extent(bb), snap='out') plot(srtm) ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-57-1.png" style="display: block; margin: auto;" /> [CIAT](http://srtm.csi.cgiar.org) --- ## Slope ```r slope <- terrain(srtm, opt='slope') plot(slope) ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-58-1.png" style="display: block; margin: auto;" /> --- ## Aspect ```r aspect <- terrain(srtm, opt='aspect') plot(aspect) ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-59-1.png" style="display: block; margin: auto;" /> --- ## Hill Shade ```r hs <- hillShade(slope, aspect, 10, 217.92) plot(hs, col= heat.colors(25)) ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-60-1.png" style="display: block; margin: auto;" /> [Azimuth Angle](https://www.esrl.noaa.gov/gmd/grad/solcalc/) --- ## Hill Shade ```r hs <- hillShade(slope, aspect, 90, 217.92) plot(hs, col= heat.colors(25)) ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-61-1.png" style="display: block; margin: auto;" /> --- ## Contour Lines ```r hs <- hillShade(slope, aspect, 45, 217.92) plot(hs, col= heat.colors(25)) omahaContour <- rasterToContour(srtm, maxpixels = 500) plot(omahaContour, add = TRUE) ``` <img src="rgis_tutorial_website_files/figure-html/unnamed-chunk-62-1.png" style="display: block; margin: auto;" /> --- ## Shiny Web Applications - Lightweight framework for web applications - Use R engine for analysis - R Interface to producing HTML, Javascript for User Interface - Customizable HTML, CSS, and Javascript ```r install.packages('shiny') install.packages('shinythemes') install.packages('osrm') ``` [Shiny](https://shiny.rstudio.com/) --- ## Open Source Routing Machine ```r library(osrm) driveTime <- osrmIsochrone(loc = c(-95.9980, 41.2524)) driveTime ``` ``` class : SpatialPolygonsDataFrame features : 6 extent : -96.92043, -95.16736, 40.522, 41.94424 (xmin, xmax, ymin, ymax) coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 variables : 4 names : id, min, max, center min values : 1, 0, 10, 5 max values : 6, 50, 60, 55 ``` [osrm](https://cran.r-project.org/web/packages/osrm/index.html), [Project OSRM](http://project-osrm.org/) --- ## Self Service Isochrone - User Input - Coordinates - Drive Times - Background map - Output - map of drive times - ability to download shapefiles --- ## Shiny Template ```r ui <- fluidPage() server <- function(input, output) { } shinyApp(ui, server) ``` - reactive programming environment - `input` is provided by the end user - `output` is processed by the server based on `input` and rendered to the end user - single file application --- class: inverse, middle, center ## Questions?