---
title: "PostGIS:"
author: "<i>Oh, the Places You'll Index!</i><br><br><br>"
date: "<small><a href='https://jacksonvoelkel.com'>Jackson Voelkel</a><br><br>Kaiser Permanente&nbsp;&nbsp;&nbsp;|&nbsp;&nbsp;&nbsp;Portland State University<br><br></small>"
output:
  revealjs::revealjs_presentation:
    theme: "simple"
    transition: slide
    self_contained: false
    reveal_plugins: ["zoom"]
editor_options: 
  chunk_output_type: console
---
 
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, results = 'none', warning = FALSE, error = FALSE, message = FALSE)
options(tigris_use_cache = TRUE)
```

## Overview

1. Who am I
2. What is GIS
3. What are Spatial Analytics
4. What is PostGIS, and where does it fit in?
5. PostGIS Examples

# Who am I?

<!-- ## Who am I? -->

<!-- <br> -->

<!-- <img src='./img/ex_crisis.png'><br> -->

## Who am I?

<br>

<p class='fragment' style='align-text:left'>Jackson Voelkel</p>
<br>
<p class='fragment' style='align-text:left'><a href='https://jacksonvoelkel.com' target='_blank'>jacksonvoelkel.com</a></p>
<br><br>
<p class='fragment' style='font-size:13pt;align-text:left'><a href='https://latlon.io' target='_blank'>I've also been saying for years I'll start a blog...</a></p>

## Health Data Analyst, Kaiser Permanente (National GIS team)

<img src = 'https://jacksonvoelkel.com/geog597/GEOG579_W2019__-_1_-_Introduction_files/img/gisatkp.jpg' ></img>

- Spatial statistics and custom analysis
- Workflow automation
- Database design/compilation
- R, Python, SQL

## Adjunct Professor of Geography, Portland State University

<img src = './img/psu.png'>

- GEOG 597: Advanced Spatial Quantitative Analysis in R
- GEOG 575: Spatial Database Design and Compilation
- GEOG 4/588: GIS 1: Introduction to GIS

## Principal / Lead Analyst, Lone Fir Labs (Consulting)

<img src = 'https://jacksonvoelkel.com/geog597/GEOG579_W2019__-_1_-_Introduction_files/img/lfl_logo.jpg' ></img>

- Same as above, but far more specialized/specific tasks
- Remote Sensing
- Machine Learning
- Web map application development


<!-- ## Health Data Analyst, Kaiser Permanente (National GIS team) -->

<!-- <img src = 'https://jacksonvoelkel.com/geog597/GEOG579_W2019__-_1_-_Introduction_files/img/gisatkp.jpg' ></img> -->

<!-- - Spatial statistics and custom analysis -->
<!-- - Workflow automation -->
<!-- - Database design/compilation -->
<!-- - R, Python, SQL -->

<!-- ## Principal / Lead Analyst, Lone Fir Labs (Consulting) -->

<!-- <img src = 'https://jacksonvoelkel.com/geog597/GEOG579_W2019__-_1_-_Introduction_files/img/lfl_logo.jpg' ></img> -->

<!-- - Same as above, but far more specialized/specific tasks -->
<!-- - Remote Sensing -->
<!-- - Machine Learning -->
<!-- - Web map application development -->

<!-- ## How did I get into that? -->

<!-- <br> -->

<!-- - Art School Dropout (I'm still into photography, but mostly images from satellites) -->
<!-- - Bachelor of Science: Geography -->
<!-- - Master of Urban Studies -->
<!-- <p class='fragment' style='test-align:center'><i>Urban Geocomputation:  -->
<!-- Two Studies on Urban Form and its Role in Altering Climate</i> <a href = "https://pdxscholar.library.pdx.edu/open_access_etds/4350/" target="_blank">[link]</a></p> -->

<!-- ## Previous Research -->

<!-- <br><br><br> -->
<!-- <small>Just to name a few broad projects...</small> -->

<!-- ## Trees, LiDAR -->

<!-- <br> -->

<!-- <img height='100%' src='https://jacksonvoelkel.com/geog597/GEOG579_W2019__-_1_-_Introduction_files/img/psu_lib.png'></img> -->

<!-- ## -->

<!-- <br> -->

<!-- <img src='https://jacksonvoelkel.com/geog597/GEOG579_W2019__-_1_-_Introduction_files/img/dt_canopy.jpg'></img> -->

<!-- ## Urban Heat Island Effect -->

<!-- <div class='stretch'> -->
<!-- <iframe class="stretch" style="height:100%;width:100%; z-index:999999999;" data-src="https://climatecope.research.pdx.edu/dc/uhi/wm_pm">[https://climatecope.research.pdx.edu/dc/uhi/wm_pm]</iframe> -->
<!-- </div> -->

<!-- <a href='https://www.mdpi.com/2225-1154/7/1/5' target='_blank'><small><i>Integrating Satellite and Ground Measurements for Predicting Locations of Extreme Urban Heat</i>, Shandas, Voelkel, Williams, and Hoffman</small></a> -->

<!-- ##  -->

<!-- <img src='./img/nyt.png'> -->
<!-- <small style='font-size:13pt'>Though they didn't even mention my name... <a target='_blank' href='https://www.nytimes.com/interactive/2019/08/09/climate/city-heat-islands.html'>New York Times, August 2019</a></small> -->

<!-- ##  -->

<!-- <img src='./img/wapo.png'> -->
<!-- <small style='font-size:13pt'>They (sort of) remembered me... <a target='_blank' href='https://www.washingtonpost.com/news/capital-weather-gang/wp/2018/07/11/noaa-heat-island-campaign-to-map-what-areas-in-d-c-and-baltimore-swelter-the-most/'>Washington Post, August 2018</a></small> -->

<!-- ## Healthcare -->

<!-- <br> -->

<!-- - Wildfire home damage and member locations -->
<!-- - Health Indices -->
<!-- - Social Determinants of Health -->
<!-- - Data, data, data -->

<!-- ## -->

<!-- <br><br> -->

<!-- But enough about that... -->

# What is GIS?

## That's a hard question to answer

<br>

- It depends on who you ask
- It depends on the audience's interests

## Common Responses:

<br>

- "Digital Cartography"
- "Math with maps"
- "Spatial data management"
- "Geodatabases"

<p class='fragment'>... or you can just watch a video.</p>

##

<iframe width="100%" height="100%" class='stretch' src="https://www.youtube.com/embed/LHDCRjAxpI0" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
<small style='font-size:13pt'>Esri is the main GIS software company around today.</small>

## 

<iframe width="100%" height="100%" class='stretch' src="https://www.youtube.com/embed/cZia3ShzTWM?start=30" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

<small>Here's William Shatner talking about one of the first GIS tools!</small>

##

<br>

<img src="./img/shat.jpg">

## Why is GIS Important?

<br>

<h2>Everything happens somewhere.</h2>

<br>

<p class='fragment'>Geographic space is often the lowest common denomenator for compiling and analyzing disparate data.</p>

# Spatial Analysis

> Spatial Analysis is concerned with investigating the <b>patterns</b> that arise as a result of <b>processes</b> that may be operating in space. Techniques and methods to enable the representation, description, measurement, comparison and generation of spatial patterns are central to the study of geographic information analysis.

<small>O'Sullivan and Unwin, page 4</small>

## Four Types of Spatial Analysis

<br>

1. Spatial data manipulation
2. Spatial data analysis
3. Spatial statistical analysis
4. Spatial modeling

## Spatial data manipulation

<br>

- Typically done via GIS
- E.g., merging, overlay, dissolve.

## Spatial data analysis

<br>

First steps in analysis: <b>exploratory</b> and <b>descriptive</b>

## Spatial statistical analysis

<br>

Hypothesis testing

## Spatial Modeling

<br>

- Forecasts and predictions
- E.g., climate change modeling

# Spatial Data Analysis

## {data-transition='none'}

<blockquote style='font-size:20pt'>Spatial data analysis is concerned with that branch of data analysis where the
geographical referencing of objects contains important information. In many areas of data
collection, especially in some areas of experimental science, the indexes that distinguish
different cases can be exchanged without any loss of information. All the information
relevant to understanding the variation in the data set is contained in the observations and
no relevant information is contained in the indexing. In the case of spatial data the
indexing (by location and time) may contain crucial information. A definition of spatial
analysis (of which spatial data analysis is one element) is that it represents a collection of
techniques and models that explicitly use the spatial referencing of each data case.</blockquote>

<small>Goodchild, Michael, and Haining. 2004. "GIS and spatial data analysis: converging perspectives". Papers in Regional Science 83:363-385</small>


## {data-transition='none'}

<blockquote style='font-size:20pt'>Spatial data analysis is concerned with that branch of data analysis where the
geographical referencing of objects contains important information. In many areas of data
collection, especially in some areas of experimental science, the indexes that distinguish
different cases can be exchanged without any loss of information. All the information
relevant to understanding the variation in the data set is contained in the observations and
no relevant information is contained in the indexing. In the case of spatial data the
indexing (by location and time) may contain crucial information. A definition of spatial
analysis (of which spatial data analysis is one element) is that it represents <b><span style='color:red'>a collection of
techniques and models that explicitly use the spatial referencing of each data case</span></b>.</blockquote>

<small>Goodchild, Michael, and Haining. 2004. "GIS and spatial data analysis: converging perspectives". Papers in Regional Science 83:363-385</small>

## How is it different than GIS?

<br>

Different perspective on spatial data

<br>

<p class='fragment'>GIS is built around the entity-attribute model. Spatial analysis uses these data, and <b>sees data as patterns</b> that are the outcomes of processes. This is similar... but different than spatial data manipulation.</p>

## Proprietary Tools

<br>

- ArcGIS
- ENVI
- MicroStation

## FOSS Tools

<br>

- R
- Python
- QGIS
- GRASS (still around, and amazing!)

# Enterprise Spatial Analysis

## Arcgis

<br>

- They 'own' GIS
- Desktop/Server software
- Python API
- Some DB development (though nothing like PostgreSQL)

## Example 1: Enriching data

We wish to filter a specific cohort of people by distance to rivers, perform and analysis, and push the results back to our DB.

1. Create Oracle Query to grab cohort of interest (<b>not filtered by distance to rivers!</b>)
2. Run query and download as .csv
3. Open .csv in ArcGIS
4. Create points from XY columns
5. Run analysis (buffer, spatial join, etc)
6. Export results to csv 
7. Upload to Oracle

## Example 2: Joining data

We wish to join the newest Census Blocks (~12 million) to every point we have in a database (lets imagine ~12 million patients, trees, or aardvarks).

1. Select all XY from Oracle
2. Download a CSV
3. Import into ArcGIS
4. Join
5. <b>Wait many hours</b>.
6. Export results as a CSV
7. Upload to Oracle.

## Aside: Oracle Spatial

<br>

Yes, I could do some of this in Oracle Spatial, but I'm not made of money <span class='fragment'>... anyways, I think PostGIS is better!</span>

## Limits

<br>

- Time: this takes <i>forever</i>!
- Reproducability: too much hands-on work
- Dependencies: Between Oracle and ArcGIS there are enough bugs to drive you mad... imagine when they update!
- Versioning: Floating files of different vintages

# Where Does PostgreSQL Fit In?

## Everywhere!

## Storage

<br>

Well, obviously!

## (Spatial) Analysis

<br>

PostGIS!

# PostGIS: Spatial Extention

##

<br>

As easy as...

<br>

```
$ sudo apt install postgresql-10-postgis-2.4
```
<small style='font-size:13'>Or the equivalent for your inferior (just kidding) system/version</small>

and <br>

```
CREATE EXTENSION POSTGIS;
```

## Adds in Massive Amounts of Power

<br>

- Data types
- Analytical functions
- So big, it has more extensions for processes like volumetrics and street network analysis.

## Geometry

Spatial data works in a standard entity-attribute format, with various geometry types.

<img src='./img/ent_at.png'>

## Geometry

The Open Geographic Consortium sets standards, which PostGIS follows.

<img src='./img/ogc_std.png'>

<small style='font-size:13pnt'><a href='https://www.opengeospatial.org/standards/sfa'>OGC Simple Feature Access - Part 1: Common Architecture</a></small>

## Geometry

```{r echo=F}
source('/var/www/html/presentations/.connections.R')
```

Let's look at a very simple data set: fire stations in the Portland Metro Area:

```{r echo=F}
if(!file.exists("/var/www/html/presentations/fire_station_map.html")){
  require(leaflet)
  require(sf)
  m <- leaflet() %>% addProviderTiles("Stamen.Toner") %>% addMarkers(data=st_read(conn,query="select * from fire_sta;") %>% st_cast("POINT") %>% st_transform(4326))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/fire_station_map.html",selfcontained = T)
}
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./fire_station_map.html'></iframe>
</div>

##

```SQL
SELECT * FROM fire_sta;
```

<pre class='fragment' style='overflow:auto'><span class="inner-pre" style="font-size:13px"> id  |                                 geom                                 | station |               address               |      city       |               district               
-----+----------------------------------------------------------------------+---------+-------------------------------------+-----------------+--------------------------------------
   1 | 01040000201E6900000100000001010000006848EEED96711F410D8911649E305341 | 50      | 12617 SW WALNUT ST                  | TIGARD          | TUALATIN VALLEY FIRE & RESCUE
   2 | 01040000201E6900000100000001010000003057F1DB66D61F41779684954A2D5341 | 211     | 4555 JEAN RD                        | LAKE OSWEGO     | LAKE OSWEGO F&R&LIFE SAFETY
   3 | 01040000201E6900000100000001010000000B36E71F0FDB1D412D5B860991315341 | 11      | 102 E MAIN ST                       | GASTON          | GASTON RFPD
   4 | 01040000201E6900000100000001010000001E8D7697A0D61F41D67BD6BA5B415341 | 26      | 5247 N LOMBARD ST                   | PORTLAND        | PORTLAND FIRE & RESCUE
   5 | 01040000201E6900000100000001010000006F3DEF04B1A51F4197EB934FA6265341 | 56      | 8455 SW ELLIGSEN RD                 | WILSONVILLE     | TUALATIN VALLEY FIRE & RESCUE
   6 | 01040000201E6900000100000001010000008FE7EC0CDF06204146A2177BF1195341 | 65      | 26815 HIGHWAY 170                   |                 | CANBY RFPD #62
   7 | 01040000201E6900000100000001010000002FC149D18B392241CBD9801B12245341 | 253     | 87600 E GOVERNMENT CAMP LOOP        | GOVERNMENT CAMP | HOODLAND RFPD #74
   8 | 01040000201E690000010000000101000000B046D46B25882041D7EC8E183C3C5341 | 74      | 1520 NE 192ND AVE                   | GRESHAM         | GRESHAM FIRE & EMERG SRVCS
   9 | 01040000201E690000010000000101000000C668D7930D901F41350FEBB7F2335341 | 53      | 8480 SW SCHOLLS FERRY RD            | BEAVERTON       | TUALATIN VALLEY FIRE & RESCUE
  10 | 01040000201E69000001000000010100000082EDAD44C9E91F417153C4A654345341 | 18      | 8720 SW 30TH AVE                    | PORTLAND        | PORTLAND FIRE & RESCUE
</span></pre>

##

```SQL
SELECT  id
        , station
        , address
        , city
        , st_astext(geom) geom 
FROM fire_sta;
```

<pre class='fragment' style='overflow:auto;'><span class="inner-pre" style="font-size:13px">
 id | station |           address            |      city       |                     geom                      
----+---------+------------------------------+-----------------+-----------------------------------------------
  1 | 50      | 12617 SW WALNUT ST           | TIGARD          | MULTIPOINT(515173.732354289 5030521.56357027)
  2 | 211     | 4555 JEAN RD                 | LAKE OSWEGO     | MULTIPOINT(521625.714787829 5027114.33621751)
  3 | 11      | 102 E MAIN ST                | GASTON          | MULTIPOINT(489155.781155438 5031492.14882545)
  4 | 26      | 5247 N LOMBARD ST            | PORTLAND        | MULTIPOINT(521640.147913175 5047662.91934105)
  5 | 56      | 8455 SW ELLIGSEN RD          | WILSONVILLE     | MULTIPOINT(518508.254818878 5020313.24340334)
  6 | 65      | 26815 HIGHWAY 170            |                 | MULTIPOINT(525167.525244938 5007301.9233175)
  7 | 253     | 87600 E GOVERNMENT CAMP LOOP | GOVERNMENT CAMP | MULTIPOINT(597189.908765828 5017672.42973943)
  8 | 74      | 1520 NE 192ND AVE            | GRESHAM         | MULTIPOINT(541714.710603913 5042416.38372346)
  9 | 53      | 8480 SW SCHOLLS FERRY RD     | BEAVERTON       | MULTIPOINT(517123.394376409 5033930.87372189)
 10 | 18      | 8720 SW 30TH AVE             | PORTLAND        | MULTIPOINT(522866.317069732 5034322.60573278)
</span></pre>

<small class='fragment' >Notice the projection of the data gives hard-to-read coordinates!</small>

##

```SQL
SELECT  id
        , station
        , address
        , city
        , (ST_ASTEXT(ST_TRANSFORM(geom,4326))) geom 
FROM fire_sta;
```

<pre class='fragment' style='overflow:auto'><span class="inner-pre" style="font-size:13px">
 id | station |           address            |      city       |                      geom                      
----+---------+------------------------------+-----------------+------------------------------------------------
  1 | 50      | 12617 SW WALNUT ST           | TIGARD          | MULTIPOINT(-122.806027592858 45.428051325105)
  2 | 211     | 4555 JEAN RD                 | LAKE OSWEGO     | MULTIPOINT(-122.723699509766 45.3972124570865)
  3 | 11      | 102 E MAIN ST                | GASTON          | MULTIPOINT(-123.138647942155 45.4368683505389)
  4 | 26      | 5247 N LOMBARD ST            | PORTLAND        | MULTIPOINT(-122.722608693335 45.5821707668332)
  5 | 56      | 8455 SW ELLIGSEN RD          | WILSONVILLE     | MULTIPOINT(-122.763784232417 45.3360833873664)
  6 | 65      | 26815 HIGHWAY 170            |                 | MULTIPOINT(-122.679455050874 45.2187565962888)
  7 | 253     | 87600 E GOVERNMENT CAMP LOOP | GOVERNMENT CAMP | MULTIPOINT(-121.760251085664 45.3058277709189)
  8 | 74      | 1520 NE 192ND AVE            | GRESHAM         | MULTIPOINT(-122.46574190944 45.5340339071676)
  9 | 53      | 8480 SW SCHOLLS FERRY RD     | BEAVERTON       | MULTIPOINT(-122.780985680916 45.4586940268239)
 10 | 18      | 8720 SW 30TH AVE             | PORTLAND        | MULTIPOINT(-122.707514263889 45.462055584884)
</span></pre>


## Polygons

Here we have the Neighborhoods in the Portland Metro Area:

```{r echo=F}
if(!file.exists("/var/www/html/presentations/fire_station_map.html")){
  require(leaflet)
  require(sf)
  nbo <- st_read(conn,query="select * from nbo_hood;") %>% st_transform(4326) %>% st_cast("POLYGON")
  col_func <- colorFactor(sample(rainbow(26),26), domain = nbo$name)
  nbo$nbo_cols <- col_func(nbo$name)
  m <- leaflet() %>% 
    addProviderTiles("Stamen.Toner") %>% 
    addPolygons(
      data=nbo, 
      weight = 1, 
      color = 'black', 
      fillColor = nbo$nbo_cols, 
      popup = paste0("<b>Name</b>: ",nbo$name))
  m2 <- leaflet() %>% 
    addProviderTiles("Stamen.Toner") %>% 
    addPolygons(
      data=nbo, 
      weight = 1, 
      color = 'black', 
      fillOpacity = 0.7, 
      fillColor = nbo$nbo_cols, 
      popup = paste0("<b>Name</b>: ",nbo$name))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/nbo_map.html",selfcontained = T)
  htmlwidgets::saveWidget(m2,"/var/www/html/presentations/nbo2_map.html",selfcontained = T)
}
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./nbo_map.html'></iframe>
</div>

##

```SQL
SELECT  nbo_hood_id
        , name
        , ST_ASTEXT(ST_TRANSFORM(geom,4326)) 
FROM nbo_hood;
```

<pre class='fragment' style='overflow:auto'><span class="inner-pre" style="font-size:13px"> nbo_hood_id |          name           |                                                                                                                                                                                           
-------------+-------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
           8 | LLOYD DISTRICT          | MULTIPOLYGON(((-122.650480796952 45.5333589655569,-122.650510431462 45.5314726993377,-122.654567304327 45.531497367206,-122.654583838372 45.5293616414892,-122.653311176239 45.52876729323
         211 | GRESHAM - NORTH GRESHAM | MULTIPOLYGON(((-122.460068956153 45.5591750869026,-122.460071836589 45.5581701759594,-122.460071474271 45.5581033059463,-122.459943494087 45.5581034443263,-122.459921124039 45.5581034656
          96 | PARKROSE HEIGHTS        | MULTIPOLYGON(((-122.537634551912 45.5337757045262,-122.537636878715 45.5335736333397,-122.537824367438 45.5335715422847,-122.545188064426 45.533490336749,-122.547000762669 45.53347027081
          28 | NORTHFIELD              | MULTIPOLYGON(((-122.497756840315 45.6083488130372,-122.495572604578 45.6083495710326,-122.495576108117 45.6101271754758,-122.493112204749 45.6101643722772,-122.49309425842 45.60917751321
         287 | NORTH SALMON CREEK      | MULTIPOLYGON(((-122.655676626356 45.7289193353581,-122.655477162054 45.7282245128809,-122.654758732645 45.7257234828009,-122.654380547849 45.7247388828086,-122.654560664795 45.7235452739
         154 | CASCADE HIGHLANDS       | MULTIPOLYGON(((-122.503939226142 45.6102082210104,-122.503923730676 45.6072320292005,-122.503908694318 45.6043453309102,-122.504661932669 45.6043443269167,-122.504804951482 45.6043411063
         162 | FISHERS LANDING EAST    | MULTIPOLYGON(((-122.495576108117 45.6101271754758,-122.495572604578 45.6083495710326,-122.497756840315 45.6083488130372,-122.497753173955 45.6066773378049,-122.492101178464 45.6066973944
         350 | CHILDS                  | MULTIPOLYGON(((-122.712797021563 45.3873736667375,-122.71278317843 45.3862934000188,-122.712778032241 45.3859530356456,-122.712773822213 45.3856748658333,-122.712769287943 45.38535516713
         147 | FOURTH PLAIN VILLAGE    | MULTIPOLYGON(((-122.639197560032 45.6458942572195,-122.638937368859 45.6458664162537,-122.638389870996 45.6463356847865,-122.63783388227 45.646320363385,-122.637376439233 45.646147518485
         316 | WEST HAZEL DELL         | MULTIPOLYGON(((-122.714306015297 45.6971360464785,-122.713822245802 45.6969213359745,-122.706154919092 45.6970914144172,-122.703045332396 45.6971602242795,-122.702985164086 45.6971605011
</span></pre>

<small style='font-size:12pt'>I've cut off the entire polygon description.</small>

# PostGIS Analysis

##

<br>

Ok, so we can store geographic information... let's put it to work!

## Spatial Functions

<small>Standard GIS operates with "tools"/"widgets":</small>

<div class="stretch" style="height:60%">
<div style='float:left;display:inline;width:50%;'>
<small>ArcGIS Toolbox</small><br>
<img src='./img/arctoolbox.png'>
</div>
<div style='float:left;display:inline;width:50%'>
<small>QGIS Toolbox</small><br>
<img style='height:40%' src='./img/qgistoolbox.png'>
</div>
</div>

## 

In PostGIS, we have functions (lots of them):

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='https://postgis.net/docs/PostGIS_Special_Functions_Index.html'></iframe>
</div>

## Post GIS Functions

<br>

- PostGIS has a massive amount of high-performance tools
- Between PostGIS and R/Python, you can do anything
- Analysis can take place where the data lives - no unnecessary transfers!

##

<br>

<p>The next few examples are <i>very basic</i>! PostGIS is capable of highly advanced analytics, which we wont get into here.</p>

# Buffering

## Buffering

<small>Here we have a cohort of schools:</small>

```{r echo=F}
if(!file.exists("/var/www/html/presentations/school1_map.html")){
  school <- sf::st_read(conn, query = "SELECT * FROM schools WHERE grade = '9-12';")
  require(leaflet)
  require(sf)
  m <- leaflet() %>% addProviderTiles("Stamen.TonerLite") %>% addMarkers(data=school %>% st_cast("POINT") %>% st_transform(4326))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/school1_map.html",selfcontained = T)
}
```

```SQL
SELECT * FROM schools WHERE grade = '9-12';
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./school1_map.html'></iframe>
</div>

## Buffering

<small>Here we have a buffer for each (converting to different coordinate system):</small>

```{r echo=F}
if(!file.exists("/var/www/html/presentations/school2_map.html")){
  school_buf <- sf::st_read(conn, query = "SELECT st_buffer(st_transform(geom,26910),2000) FROM schools WHERE grade = '9-12';") %>% st_cast("POLYGON")
  require(leaflet)
  require(sf)
  m <- leaflet() %>% addProviderTiles("Stamen.TonerLite") %>% addPolygons(data=school_buf %>% st_transform(4326), weight = 1, color = 'black')
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/school2_map.html",selfcontained = T)
}
```

```SQL
SELECT ST_BUFFER(ST_TRANSFORM(geom,26910),2000) 
FROM schools WHERE grade = '9-12';
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./school2_map.html'></iframe>
</div>

## Buffer + Union

<small>Here we have a single 'unioned' buffer:</small>

```{r echo=F}
if(!file.exists("/var/www/html/presentations/school3_map.html")){
  require(leaflet)
  school_buf_dis <- sf::st_read(conn, query = "SELECT st_union(st_buffer(st_transform(geom,26910),2000)) FROM schools WHERE grade = '9-12';") %>% st_cast("POLYGON")
  require(sf)
  m <- leaflet() %>% addProviderTiles("Stamen.TonerLite") %>% addPolygons(data=school_buf_dis %>% st_transform(4326), weight = 2, color = 'red',fillColor = 'red')
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/school3_map.html",selfcontained = T)
}
```

```SQL
SELECT ST_UNION(ST_BUFFER(ST_TRANSFORM(geom,26910),2000)) 
FROM schools WHERE grade = '9-12';
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./school3_map.html'></iframe>
</div>

# Spatial Join

## Spatial Join

<br>

One of the most common functions in GIS is a "spatial join", wherein tables are joined based on geographic relation. 

- In most desktop software this can take hours on large data.
- In PostGIS, it is blazingly fast.

<p class='fragment'>Let's join the neighborhood name onto each school in the region!</p>

## Spatial Join

<small>As a reminder, here are the schools and neighborhoods:</small>

<div class='stretch'>
<iframe class='stretch' height='100%' width='49%' src='./nbo_map.html'></iframe>
<iframe class='stretch' style='float:left;height:100%;width:49%' src='./school1_map.html'></iframe>
</div>

## Spatial Join

```{r echo=F, eval=F, warning=FALSE, message=FALSE, error=FALSE}
if(!file.exists("/var/www/html/presentations/join_map.html")){
  schl_nbo <- sf::st_read(conn, query = "SELECT s.*, n.name nbo_name FROM schools s JOIN nbo_hood n ON ST_Intersects(s.geom,st_transform(n.geom,4326)) where s.grade = '9-12';") %>% st_cast("POINT")
  require(leaflet)
  require(sf)
  m <- leaflet() %>% 
    addProviderTiles("Stamen.TonerLite") %>% 
    addCircleMarkers(
      data=schl_nbo, 
      weight = 1, 
      color = 'black', 
      fillOpacity = 0.7, 
      fillColor = st_join(schl_nbo,nbo)$nbo_cols , 
      popup = paste0("<b>School</b>: ",schl_nbo$name,"<br><b>NBO</b>: ", schl_nbo$nbo_name))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/join_map.html",selfcontained = T)
}
```

```SQL
SELECT s.*, n.name nbo_name 
FROM schools s 
JOIN nbo_hood n 
  ON ST_Intersects(s.geom,st_transform(n.geom,4326)) 
WHERE s.grade = '9-12'; 
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='49%' src='./nbo2_map.html'></iframe>
<iframe class='stretch' style='float:left;height:100%;width:49%' src='./join_map.html'></iframe>
</div>

## Spatial Join: Point-in-Polygon

# Intersections

## Intersections

<br>

Let's get all of the Census Block Groups (CBGs) within Downtown Portland. There are many options for this:

- "Clip" CBGs to neighborhood
- Get all CBGs with centroids in the neighborhood
- Get all CBGs that touch or 'intersect' the neighborhood

## Intersections: 'Clipping'

```{r echo = F, warning=FALSE, error=FALSE}
if(!file.exists("/var/www/html/presentations/int_clip.html")){
  int_clip <- sf::st_read(conn, query = "SELECT c.fips, c.pop10, ST_Intersection(n.geom, c.geom) geom FROM cbg c, nbo_hood n WHERE n.name = 'DOWNTOWN' AND ST_Intersects(n.geom,c.geom);")
  require(leaflet);require(sf)
  m <- leaflet() %>% 
    addProviderTiles("Stamen.TonerLite") %>% 
    addPolygons(
      data= int_clip %>% st_transform(4326), 
      weight = 5,
      color = 'white',
      opacity=1,
      fillColor = 'blue',
      fillOpacity = 0.3,
      popup = paste0("<b>Population (2010)</b>: ", int_clip$pop10))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/int_clip.html",selfcontained = T)
}
```

```SQL
SELECT  c.fips
        , c.pop10
        , ST_Intersection(n.geom, c.geom) geom
FROM  cbg c, 
      nbo_hood n 
WHERE n.name = 'DOWNTOWN'
  AND ST_Intersects(n.geom,c.geom);
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./int_clip.html'></iframe>
</div>

## Intersections: 'Touching'

```{r echo = F, warning=FALSE, error=FALSE}
if(!file.exists("/var/www/html/presentations/int_touch.html")){
  int_touch <- sf::st_read(conn, query = "SELECT c.* FROM  cbg c, nbo_hood n WHERE ST_INTERSECTS(c.geom, n.geom) AND n.name = 'DOWNTOWN';")
  require(leaflet);require(sf)
  m <- leaflet() %>% 
    addProviderTiles("Stamen.TonerLite") %>% 
    addPolygons(
      data= int_touch %>% st_transform(4326), 
      weight = 4,
      color = 'white',
      opacity=1,
      fillColor = 'blue',
      fillOpacity = 0.3,
      popup = paste0("<b>Population (2010)</b>: ", int_clip$pop10))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/int_touch.html",selfcontained = T)
}
```

```SQL
SELECT c.*
FROM  cbg c, 
      nbo_hood n 
WHERE ST_INTERSECTS(c.geom, n.geom) 
  AND n.name = 'DOWNTOWN';
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./int_touch.html'></iframe>
</div>

## Intersections: 'Centroids'

```{r echo = F, warning=FALSE, error=FALSE}
if(!file.exists("/var/www/html/presentations/int_cent.html")){
  int_cent <- sf::st_read(conn, query = "SELECT c.* FROM  cbg c, nbo_hood n WHERE ST_INTERSECTS(ST_CENTROID(c.geom), n.geom) AND n.name = 'DOWNTOWN';")
  require(leaflet);require(sf)
  m <- leaflet() %>% 
    addProviderTiles("Stamen.TonerLite") %>% 
    addPolygons(
      data= int_cent %>% st_transform(4326), 
      weight = 4,
      color = 'white',
      opacity=1,
      fillColor = 'blue',
      fillOpacity = 0.3,
      popup = paste0("<b>Population (2010)</b>: ", int_clip$pop10))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/int_cent.html",selfcontained = T)
}
```

```SQL
SELECT c.*
FROM  cbg c, 
      nbo_hood n 
WHERE ST_INTERSECTS(ST_CENTROID(c.geom), n.geom) 
  AND n.name = 'DOWNTOWN';
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./int_cent.html'></iframe>
</div>

## Intersections: Density

```{r echo = F, warning=FALSE, error=FALSE}
if(!file.exists("/var/www/html/presentations/pop_dens.html")){
  pop_dens <- sf::st_read(conn, query = "SELECT fips, pop10/(ST_AREA(c.geom)/1000) pop_density, c.geom FROM  cbg c, nbo_hood n WHERE ST_INTERSECTS(ST_CENTROID(c.geom), n.geom) AND n.name = 'DOWNTOWN';")
  require(leaflet)
  require(sf)
  dens_colors <- colorQuantile("RdBu", domain = pop_dens$pop_density, n=7,reverse = T)
  m <- leaflet() %>% 
    addProviderTiles("Stamen.TonerLite") %>% 
    addPolygons(
      data=pop_dens %>% st_transform(4326), 
      weight = 1,
      color = 'black',
      fillOpacity = 0.7,
      fillColor = dens_colors(pop_dens$pop_density),
      popup = paste0("<b>People/km^2</b>: ", round(pop_dens$pop_density,4)))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/pop_dens.html",selfcontained = T)
}
```

<small>While we're at it, lets also use ST_AREA to create a population density map!</small>

```SQL
SELECT  fips
        , pop10/(ST_AREA(c.geom)/1000) pop_density
        , c.geom
FROM  cbg c, 
      nbo_hood n 
WHERE ST_INTERSECTS(ST_CENTROID(c.geom), n.geom) 
  AND n.name = 'DOWNTOWN';
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./pop_dens.html'></iframe>
</div>

## Intersections

<br>

These aren't the only options! There are far more options such as:

- Completely within
- Within a distance
- Share (or don't share) edges

# Length

## Trails

Let's find the top 10 longest trails in the Portland Metro Area!

```{r echo=FALSE, warning=FALSE, error=FALSE, message=FALSE}
if(!file.exists("/var/www/html/presentations/all_trails.html")){
  require(sf);require(leaflet)
  trails <- st_read(conn,query="SELECT * FROM trails WHERE ST_LENGTH(geom) > 1000;",geometry_column='geom')
  m <- leaflet() %>% addProviderTiles("Stamen.Terrain") %>% addPolylines(data=trails %>% st_transform(4326),color = "red", weight = 3, opacity = 1)
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/all_trails.html")
}
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./all_trails.html'></iframe>
</div>

## Trails

```{r echo=FALSE, warning=FALSE, error=FALSE, message=FALSE}
if(!file.exists("/var/www/html/presentations/longest_trails.html")){
  require(sf);require(leaflet)
  trails_10 <- st_read(conn,query="SELECT systemname, sum(st_length(geom))/5280 length, ST_UNION(geom) FROM trails WHERE systemname IS NOT NULL GROUP BY systemname ORDER BY length DESC LIMIT 10;")
  trail_10_col <- colorQuantile(colorRampPalette(c('yellow3','orange','darkorange','red','firebrick'))(10),domain = trails_10$length,n = 10)
  m <- leaflet() %>% 
    addProviderTiles("Stamen.Terrain") %>% 
    addPolylines(
      data=trails_10 %>% st_transform(4326),
      color = trail_10_col(trails_10$length), 
      weight = 5, 
      opacity = 1, 
      label = lapply(paste0("<b>System Name</b>: ",trails_10$systemname,"<br><b>Length</b>: ",prettyNum(round(trails_10$length,1),big.mark = ",")," miles"),htmltools::HTML))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/longest_trails.html")
}


```

```SQL
SELECT  systemname
        , sum(st_length(geom))/5280 length
        , ST_UNION(geom) -- merging them
FROM trails WHERE systemname IS NOT NULL
GROUP BY systemname ORDER BY length DESC LIMIT 10;
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./longest_trails.html'></iframe>
</div>

# Distance

## Firestations and Schools

What is the closest fire station to Riverdale High?

```SQL
SELECT  f.*,
        ST_DISTANCE(f.geom,ST_TRANSFORM(s.geom,26910)) dist
FROM  fire_sta f, 
      (SELECT * FROM schools WHERE name = 'Riverdale High') s
ORDER BY dist
LIMIT 1;
```

<pre class='fragment' style='overflow:auto'><span class="inner-pre" style="font-size:13px"> id | station |                   st_astext                   |       dist       
----+---------+-----------------------------------------------+------------------
 97 | 10      | MULTIPOINT(524763.806430033 5034486.74243309) | 1032.21244503901
</span></pre>

<small class='fragment'>Not bad, but we can do more!</small>

## Distance: Line

```{r echo = F, warning=FALSE, error=FALSE}
if(!file.exists("/var/www/html/presentations/makeline1_map.html")){
  line1 <- sf::st_read(conn, query = "WITH s AS (SELECT ST_Transform(geom,26910) geom FROM schools WHERE name = 'Riverdale High')
SELECT ST_Transform(ST_Makeline(st_centroid(f.geom), ST_Centroid(s.geom)),4326) FROM s, (SELECT  f.*, ST_DISTANCE(f.geom,s.geom) dist
FROM  fire_sta f, s ORDER BY dist LIMIT 1) f;;")
  require(leaflet)
  require(sf)
  m <- leaflet() %>% 
    addProviderTiles("Stamen.TonerLite") %>% 
    addPolylines(
      data=line1, 
      weight = 3,
      opacity = 1,
      color = 'red')
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/makeline1_map.html",selfcontained = T)
}
```

```SQL
WITH s AS 
  (SELECT ST_Transform(geom,26910) geom 
  FROM schools 
  WHERE name = 'Riverdale High')
SELECT ST_Transform(ST_Makeline(st_centroid(f.geom), ST_Centroid(s.geom)),4326)
FROM s,
(SELECT  f.*, ST_DISTANCE(f.geom,s.geom) dist
FROM  fire_sta f, s
ORDER BY dist
LIMIT 1) f;
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./makeline1_map.html'></iframe>
</div>

## Distance: All of the Lines!

<br>

We can use a `LATERAL JOIN` to calculate all nearest neighbors!

## Distance: All of the Lines!

```{r echo = F, warning=FALSE, error=FALSE}
if(!file.exists("/var/www/html/presentations/makeline2_map.html")){
  line2 <- sf::st_read(conn, query = "WITH s AS (SELECT * FROM schools WHERE grade = '9-12') SELECT ST_Transform(ST_Makeline(ST_Centroid(f.geom), ST_Centroid(ST_Transform(s.geom,26910))),4326) FROM s CROSS JOIN LATERAL (SELECT  address, geom FROM fire_sta ORDER BY ST_Transform(s.geom,26910) <#> geom LIMIT 1) f;")
  require(leaflet)
  require(sf)
  fire_stations <- st_read(conn,query="SELECT * FROM fire_sta") %>% st_cast("POINT") %>% st_transform(4326)
  high_schools <- st_read(conn,query="SELECT * FROM schools WHERE grade = '9-12'") %>% st_cast("POINT")
  m <- leaflet() %>% 
    addProviderTiles("Stamen.TonerLite") %>% 
    addPolylines(
      data=line2, 
      weight = 3,
      opacity = 1,
      color = 'black') %>% 
    addCircleMarkers(
      data = fire_stations,
      color= "black",
      weight = 1,
      fillColor = "red",
      fillOpacity = 0.6
    ) %>% 
    addCircleMarkers(
      data = high_schools,
      color= "black",
      weight = 1,
      fillColor = "blue",
      fillOpacity = 0.6
    ) %>% addLegend(position="topright",colors = c('red','blue'), labels = c("Fire Stations","Schools"))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/makeline2_map.html",selfcontained = T)
}
```

```SQL
WITH s AS (SELECT * FROM schools WHERE grade = '9-12')
SELECT ST_Transform(ST_Makeline(ST_Centroid(f.geom), ST_Centroid(ST_Transform(s.geom,26910))),4326)
FROM s CROSS JOIN LATERAL
  (SELECT address, geom
  FROM fire_sta ORDER BY ST_Transform(s.geom,26910) <#> geom LIMIT 1) f;
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./makeline2_map.html'></iframe>
</div>

## Distance: Filtering

<br> 

We can use distance to features to filter our cohort based on other features:

```SQL
SELECT  s.name
        , ROUND(ST_Distance(ST_Transform(s.geom,2913),t.geom)) as dist
        , s.geom
FROM schools s, trails t 
WHERE ST_Distance(ST_Transform(s.geom,2913),t.geom) < 100;
```



# Street Networks

## What??

<br>

Load in data from regional governments, vendors, or even Open Streetmaps!

<p class ='fragment'>However, we'll need one more extension....</p>

## pgRouting


<img src='https://jacksonvoelkel.com/geog575/slide_images/pgrouting.png'>

Like PostGIS, pgRouting is an extension that adds data types and functions for the handling of network analysis. Some cool functions are:

- All Pairs Shortest Path, Floyd-Warshall Algorithm
- Shortest Path Dijkstra
- Driving Distance
- K-Dijkstra, One to Many Shortest Path
- Traveling Sales Person

## Open Streetmaps (OSM)

<br>

- Crowd sourced
- Various qualities, though good in major American/European urban areas
- Free as .osm.pbf files from <a href='https://www.geofabrik.de/' target = '_blank'>GEOFabrik</a>
- Can be directly parsed into

## OSM2PO

<br>

I primarily use OSM data because of this tool! It will parse raw OSM to networks with:

- Speed Limits
- Road Types
- Trails
- Stairwells
- Ferry lines

# Routing

## Routing: A-to-B

<small>First, lets look at the roads data, which consists of Nodes and Edges. We're mostly interested in the for visualization and nodes for analysis:</small>

<pre style='overflow:auto'><span class="inner-pre" style="font-size:13px"> id | osm_id  |        osm_name         | osm_meta | osm_source_id | osm_target_id | clazz | flags | source | target |    km     | kmh |   cost    | reverse_cost |      x1      |     y1     |      x2      |     y2     |    geom    
----+---------+-------------------------+----------+---------------+---------------+-------+-------+--------+--------+-----------+-----+-----------+--------------+--------------+------------+--------------+------------+--------
  1 | 4440143 | Northeast Hodes Street  |          |      38756004 |      29936590 |    41 |     7 |      1 |     50 |  0.059913 |  40 | 0.0014978 |    0.0014978 | -122.8929526 | 45.5297334 | -122.8922158 | 45.5296458 | 0105000...
  2 | 4440143 | Northeast Hodes Street  |          |      29936590 |      27286355 |    41 |     7 |     50 |      6 | 0.0115988 |  40 |   0.00029 |      0.00029 | -122.8922158 | 45.5296458 | -122.8920669 | 45.5296452 | 0105000...
  3 | 4440143 | Northeast Hodes Street  |          |      27286355 |      29936591 |    41 |     7 |      6 |     54 | 0.0530484 |  40 | 0.0013262 |    0.0013262 | -122.8920669 | 45.5296452 | -122.8913859 | 45.5296425 | 0105000...
  4 | 4440143 | Northeast Hodes Street  |          |      29936591 |      29936595 |    41 |     7 |     54 |     58 | 0.0657458 |  40 | 0.0016436 |    0.0016436 | -122.8913859 | 45.5296425 | -122.8905419 | 45.5296391 | 0105000...
  5 | 4440143 | Northeast Hodes Street  |          |      29936595 |      27199860 |    41 |     7 |     58 |     13 | 0.0685188 |  40 |  0.001713 |     0.001713 | -122.8905419 | 45.5296391 | -122.8896623 | 45.5296356 | 0105000...
</span></pre>

<small><b<Source</b> and <b>target</b> in the edges (above) links up with the node/intersection ID's (below)!</small>

<pre style='overflow:auto'><span class="inner-pre" style="font-size:13px"> id | clazz |  osm_id  | osm_name | ref_count | restrictions | geom_vertex  
----+-------+----------+----------+-----------+--------------+---------------
  1 |     0 | 38756004 |          |         2 |              | 010100002...
 50 |     0 | 29936590 |          |         1 |              | 010100002...
  6 |     0 | 27286355 |          |         1 |              | 010100002...
 54 |     0 | 29936591 |          |         1 |              | 010100002...
 58 |     0 | 29936595 |          |         1 |              | 010100002...
</span></pre>

## Routing: A-to-B

For this first route, let's route between these targets: 45667 and 34045

```{r echo=F,warning=F,message=F,error=F}
if(!file.exists("/var/www/html/presentations/ab_points.html")){
  require(sf);require(leaflet)
  ab_point <- st_read(conr,query="SELECT * FROM pdx_roads_node WHERE id IN (92015,251821);")
  m <- leaflet() %>% 
    addTiles('http://{s}.sm.mapstack.stamen.com/darkmatter/{z}/{x}/{y}.png',attribution = '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors | &copy; <a href="https://carto.com/attributions">CARTO</a>') %>% 
    addCircleMarkers(
      data=ab_point,
      radius = 8,
      weight = 1,
      color = 'black',
      fillColor = c('yellow','purple'), 
      fillOpacity = 0.9,
      popup = paste0("<b>ID</b>: ",ab_point$id))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/ab_points.html")
}
```

```SQL
SELECT * FROM pdx_roads_node WHERE id IN (92015,251821);
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./ab_points.html'></iframe>
</div>

## Routing: A-to-B


```{r echo=F,warning=F,message=F,error=F}
if(!file.exists("/var/www/html/presentations/ab_route.html")){
  require(sf);require(leaflet)
  ab_route <- st_read(conr,query="SELECT p.osm_name, x.cost, p.geom geom FROM pdx_roads p join (SELECT * FROM pgr_dijkstra('SELECT id, source, target, cost, reverse_cost FROM pdx_roads',92015,251821)) x ON p.id = x.edge")
  ab_point <- st_read(conr,query="SELECT * FROM pdx_roads_node WHERE id IN (92015,251821);")
  m <- leaflet() %>% 
    addTiles('http://{s}.sm.mapstack.stamen.com/darkmatter/{z}/{x}/{y}.png',attribution = '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors | &copy; <a href="https://carto.com/attributions">CARTO</a>') %>% 
    addPolylines(
      data=ab_route,
      weight = 9,
      color = 'blue',
      opacity = 0.9,
      label = lapply(paste0("<b>Street Name</b>: ",ab_route$osm_name,"<br><b>Cost</b>: ",round(ab_route$cost,5)),htmltools::HTML)) %>% 
    addCircleMarkers(
      data=ab_point,
      radius = 8,
      weight = 1,
      color = 'black',
      fillColor = c('yellow','purple'), 
      fillOpacity = 0.9,
      popup = paste0("<b>ID</b>: ",ab_point$id))
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/ab_route.html")
}
```

```SQL
SELECT   p.osm_name, x.cost, p.geom
FROM pdx_roads p join (
    SELECT * FROM pgr_dijkstra(
        'SELECT id, source, target, cost, reverse_cost FROM pdx_roads'
        ,92015
        ,251821
    )
) x ON p.id = x.edge;
```
<pre style='height:250px;overflow-y:scroll'><span class="inner-pre" style="font-size:13px">              osm_name               |   cost    |      geom            
-------------------------------------+-----------+-------------------
 Southeast Helene Street             |  0.001856 | 0105000020E6100000...
 Southeast 58th Avenue               | 0.0027587 | 0105000020E6100000...
 Southeast Concorde Place            | 0.0011906 | 0105000020E6100000...
 Southeast Imlay Avenue              | 0.0037218 | 0105000020E6100000...
 Southeast Lois Street               | 0.0059255 | 0105000020E6100000...
 Southeast Lois Street               | 0.0007229 | 0105000020E6100000...
 Southeast Lois Street               | 0.0021492 | 0105000020E6100000...
 Southeast Lois Street               | 0.0014191 | 0105000020E6100000...
 Southeast Lois Street               | 0.0026148 | 0105000020E6100000...
 Southeast Lois Street               |  0.002863 | 0105000020E6100000...
 Southeast Lois Street               | 0.0007276 | 0105000020E6100000...
 Southeast Lois Street               | 0.0009849 | 0105000020E6100000...
 Southeast Lois Street               | 0.0024699 | 0105000020E6100000...
 Southeast Lois Street               | 0.0011073 | 0105000020E6100000...
 Southeast Lois Street               |  0.001348 | 0105000020E6100000...
 Southeast Lois Street               | 0.0025831 | 0105000020E6100000...
 Southeast Lois Street               |  0.001708 | 0105000020E6100000...
 Southeast Lois Street               | 0.0015847 | 0105000020E6100000...
 Southeast Lois Street               | 0.0074722 | 0105000020E6100000...
 Southeast Lois Street               | 0.0019735 | 0105000020E6100000...
 Southeast Lois Street               | 0.0025121 | 0105000020E6100000...
 Southwest Cornelius Pass Road       | 0.0027643 | 0105000020E6100000...
 Southwest Cornelius Pass Road       | 0.0008814 | 0105000020E6100000...
 Southwest Cornelius Pass Road       | 0.0018234 | 0105000020E6100000...
 Southeast Cornelius Pass Road       | 0.0007404 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0007444 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0028978 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0009953 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0009114 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0004179 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0004582 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0006979 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0005732 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0004901 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0027577 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0004096 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0004396 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0009592 | 0105000020E6100000...
 Southwest Baseline Road             |   0.00151 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0011425 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0028631 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0007222 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0028163 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0005794 | 0105000020E6100000...
 Southwest Baseline Road             |  0.001034 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0028915 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0010045 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0058289 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0022246 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0002601 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0012276 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0001367 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0001642 | 0105000020E6100000...
 Southwest Baseline Road             | 0.0006795 | 0105000020E6100000...
 
</pre></span>

<small>And on, and on, and on...</small>

##

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./ab_route.html'></iframe>
</div>

## Routing: A-to-B

We can even get stats about the ride:

```SQL
WITH route as (SELECT p.osm_name, x.cost, st_transform(p.geom,26910) geom FROM pdx_roads p join (SELECT * FROM pgr_dijkstra('SELECT id, source, target, cost, reverse_cost FROM pdx_roads' ,92015,251821)) x ON p.id = x.edge)

SELECT  SUM(ST_Length(geom)/1000) Total_km
        , SUM(cost)*60 Total_Min
FROM route;
```
<pre style='overflow:auto'><span class="inner-pre" style="font-size:13px"> 
     total_km     | total_min 
------------------+-----------
 32.6837767909505 | 27.968016
</pre></span>

## Routing: [A,B,C...]-to-Z

<br>

We can route many-to-many if we'd like! 

The following app shows commute flows based off of the LEHD LODES data set and OSM data.

<br>

<h3><a href='https://ims.research.pdx.edu/comnet/' target='_blank'>COMNET</a></h3>

# Routing: Driving Distance

## Driving Distance

<br><br>

What if we wanted to see everywhere we could drive within <b>x</b> minutes?

<p class='fragment'>The "Driving Distance" tool gives us the ability to definite a starting point, and identify all points within a specified cost from it.</p>

## Driving Distance

```{r echo=F,message=F,warning=F,error=F}
if(!file.exists("/var/www/html/presentations/dd.html")){
  require(sf);require(leaflet)
  dd <- st_read(conr,query="SELECT  x.agg_cost as cost, x.node, x.edge, p.geom
FROM pdx_roads_node p join (
    SELECT * FROM pgr_drivingdistance(
        'SELECT id, source, target, cost, reverse_cost FROM pdx_roads'
        ,250376 -- start point, Downtown Portland
        ,(SELECT (SELECT 10.000) / (SELECT 60.000))) -- 10min
) x ON p.id = x.node;")
  dd_point <- st_read(conr,query="SELECT * FROM pdx_roads_node WHERE id = 250376;")
  gre_red <- colorRampPalette(c('forestgreen','green','yellow','orange','red','darkred'))(50)
  dd_cols <- colorQuantile(
    gre_red,
    domain = dd$cost,n = 30)
  m <- leaflet() %>% 
    addTiles('http://{s}.sm.mapstack.stamen.com/darkmatter/{z}/{x}/{y}.png',attribution = '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors | &copy; <a href="https://carto.com/attributions">CARTO</a>') %>% 
    addCircleMarkers(
      data=dd,
      radius = 5,
      weight = 0,
      fillColor = dd_cols(dd$cost), 
      fillOpacity = 0.9,
      popup = lapply(paste0("<b>Node</b>: ",dd$node,"<br><b>Cost</b>: ",paste0(round(dd$cost,2)*60), " minutes"),htmltools::HTML)) %>% 
    addMarkers(data=dd_point, popup = "Starting Point")
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/dd.html")
}
```

```SQL
SELECT  x.agg_cost as cost, x.node, x.edge, p.geom
FROM pdx_roads_node p join (
    SELECT * FROM pgr_drivingdistance(
        'SELECT id, source, target, cost, reverse_cost FROM pdx_roads'
        ,250376 -- start point, Downtown Portland
        ,(SELECT (SELECT 10.000) / (SELECT 60.000))) -- 10min
) x ON p.id = x.node;
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./dd.html'></iframe>
</div>

## Walking Distance

<br><br>

What if we wanted to see everywhere we could <b><i>walk</i></b> within <b>x</b> minutes?

## Walking Distance

```{r echo=F,message=F,warning=F,error=F}
if(!file.exists("/var/www/html/presentations/wd.html")){
  require(sf);require(leaflet)
  wd <- st_read(conr,query="SELECT  x.agg_cost as cost, x.node, x.edge, p.geom
FROM pdx_roads_node p join (
    SELECT * FROM pgr_drivingdistance(
        'SELECT id, source, target, km/5.00 as cost, reverse_cost FROM pdx_roads --!!!
        WHERE clazz > 14' --!!! OSM has road types!
        ,250376 -- start point, Downtown Portland
        ,(SELECT (SELECT 10.000) / (SELECT 60.000))) -- 10min
) x ON p.id = x.node;")
  dd_point <- st_read(conr,query="SELECT * FROM pdx_roads_node WHERE id = 250376;")
  gre_red <- colorRampPalette(c('forestgreen','green','yellow','orange','red','darkred'))(50)
  wd_cols <- colorQuantile(
    gre_red,
    domain = wd$cost,n = 9)
  m <- leaflet() %>% 
    addTiles('http://{s}.sm.mapstack.stamen.com/darkmatter/{z}/{x}/{y}.png',attribution = '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors | &copy; <a href="https://carto.com/attributions">CARTO</a>') %>% 
    addCircleMarkers(
      data=wd,
      radius = 5,
      weight = 0,
      fillColor = wd_cols(wd$cost), 
      fillOpacity = 0.9,
      popup = lapply(paste0("<b>Node</b>: ",wd$node,"<br><b>Cost</b>: ",paste0(round(wd$cost,2)*60), " minutes"),htmltools::HTML)) %>% 
    addMarkers(data=dd_point, popup = "Starting Point")
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/wd.html")
}
```

```SQL
SELECT  x.agg_cost as cost, x.node, x.edge, p.geom
FROM pdx_roads_node p join (
    SELECT * FROM pgr_drivingdistance(
        'SELECT id, source, target, km/5.00 as cost, reverse_cost FROM pdx_roads --!!!
        WHERE clazz > 14' --!!! OSM has road types!
        ,250376 -- start point, Downtown Portland
        ,(SELECT (SELECT 10.000) / (SELECT 60.000))) -- 10min
) x ON p.id = x.node;
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./wd.html'></iframe>
</div>

## Well, we don't know nodes...

<br>

We don't know which node we live by... so we can use our 'nearest neighbor' function to find it!

## Routing: Search

Using the following, it would be easy to create an application that created driving distances based on lat/lon (such as a click on a map):

```SQL
SELECT  x.agg_cost as cost, x.node, x.edge, p.geom
FROM pdx_roads_node p join (
    SELECT * FROM pgr_drivingdistance(
        'SELECT id, source, target, cost, reverse_cost FROM pdx_roads'
        , (SELECT p.id 
          FROM pdx_roads_node p, 
            (SELECT ST_SetSRID(ST_MakePoint(-122.551,45.4668),4326) as geom) s
          ORDER BY ST_DISTANCE(s.geom,p.geom) 
          LIMIT 1)
        ,(SELECT (SELECT 6.500) / (SELECT 60.000)))
) x ON p.id = x.node;
```

## Routing: Search


```{r echo=F,message=F,warning=F,error=F}
if(!file.exists("/var/www/html/presentations/search.html")){
  require(sf);require(leaflet)
  sd <- st_read(conr,query="SELECT  x.agg_cost as cost, x.node, x.edge, p.geom
FROM pdx_roads_node p join (
    SELECT * FROM pgr_drivingdistance(
        'SELECT id, source, target, cost, reverse_cost FROM pdx_roads'
        , (SELECT p.id 
          FROM pdx_roads_node p, 
            (SELECT ST_SetSRID(ST_MakePoint(-122.551,45.4668),4326) as geom) s
          ORDER BY ST_DISTANCE(s.geom,p.geom) 
          LIMIT 1)
        ,(SELECT (SELECT 6.500) / (SELECT 60.000)))
) x ON p.id = x.node;")
  sd_point <- st_read(conr,query="SELECT * FROM pdx_roads_node WHERE id = (SELECT p.id 
          FROM pdx_roads_node p, 
            (SELECT ST_SetSRID(ST_MakePoint(-122.551,45.4668),4326) as geom) s
          ORDER BY ST_DISTANCE(s.geom,p.geom) 
          LIMIT 1);")
  gre_red <- colorRampPalette(c('forestgreen','green','yellow','orange','red','darkred'))(50)
  sd_cols <- colorQuantile(
      gre_red,
    domain = sd$cost,n = 30)
  m <- leaflet() %>% 
    addTiles('http://{s}.sm.mapstack.stamen.com/darkmatter/{z}/{x}/{y}.png',attribution = '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors | &copy; <a href="https://carto.com/attributions">CARTO</a>') %>% 
    addCircleMarkers(
      data=sd,
      radius = 5,
      weight = 0,
      fillColor = sd_cols(sd$cost), 
      fillOpacity = 0.9,
      popup = lapply(paste0("<b>Node</b>: ",sd$node,"<br><b>Cost</b>: ",paste0(round(sd$cost,2)*60), " minutes"),htmltools::HTML)) %>% 
    addMarkers(data=sd_point, popup = "<b>Lat</b>: -122.551<br><b>Lon</b>: 45.4668")
  htmlwidgets::saveWidget(m,"/var/www/html/presentations/search.html")
}
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./search.html'></iframe>
</div>

## Routing

<br>

Things I didn't create an example for:

- Route from directly in front of a home, part way down edge
- Create driving distance polygons part way down edges
- Turn nodes into isochrones

## Routing: Isochrones:

<small>I'm not going over the code to do this, but it's from <a href='https://abelvm.github.io/sql/contour/' target='_blank'>Abel Vázquez</a></small>

```{r echo=F,message=F,error=F,warning=F}
if(!file.exists("/var/www/html/presentations/cont.html")){
    require(sf);require(leaflet)
    cont <- st_read(conr,query="
      WITH del AS (
      SELECT  x.agg_cost as cost, x.node, x.edge, p.geom
              FROM pdx_roads_node p join (
          SELECT * FROM pgr_drivingdistance(
              'SELECT id, source, target, cost, reverse_cost FROM pdx_roads'
              , (SELECT p.id 
                FROM pdx_roads_node p, 
                  (SELECT ST_SetSRID(ST_MakePoint(-122.551,45.4668),4326) as geom) s
                ORDER BY ST_DISTANCE(s.geom,p.geom) 
                LIMIT 1)
              ,(SELECT (SELECT 6.500) / (SELECT 60.000)))
      ) x ON p.id = x.node
      )
      SELECT * from contour_lines(
        ARRAY(SELECT geom FROM del),
        ARRAY(SELECT cost::numeric * 60 FROM del), 
        ARRAY[1,2,3,4,5,6]
      );")
    cont_point <- st_read(conr,query="SELECT * FROM pdx_roads_node WHERE id = (SELECT p.id 
            FROM pdx_roads_node p, 
              (SELECT ST_SetSRID(ST_MakePoint(-122.551,45.4668),4326) as geom) s
            ORDER BY ST_DISTANCE(s.geom,p.geom) 
            LIMIT 1);")
    gre_red <- colorRampPalette(c('forestgreen','green','yellow','orange','red','darkred'))(50)
    cont_cols <- colorQuantile(
      gre_red,
      domain = cont$`break`,n = 9)
    m <- leaflet() %>% 
      addTiles('http://{s}.sm.mapstack.stamen.com/darkmatter/{z}/{x}/{y}.png',attribution = '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors | &copy; <a href="https://carto.com/attributions">CARTO</a>') %>% 
      addPolylines(
        data=cont,
        weight = 5,
        color = cont_cols(cont$`break`), 
        opacity = 0.9,
        label = lapply(paste0("<b>Minutes</b>: ",cont$`break`),htmltools::HTML)) %>% 
      addMarkers(data=cont_point, popup = "Starting Point")
    htmlwidgets::saveWidget(m,"/var/www/html/presentations/cont.html")
}
```

<div class='stretch'>
<iframe class='stretch' height='100%' width='100%' src='./cont.html'></iframe>
</div>

# The GIS Stack

## What Does a Spatial Database Do?

<img src='https://jacksonvoelkel.com/geog575/images/data_model1.png'>

##

<img src='https://jacksonvoelkel.com/geog575/images/DBMS.jpg'>

## Esri Stack (Proprietary)

<img src='https://jacksonvoelkel.com/geog575/images/ESRI_GIS_Stack.jpg'>

##

<iframe class="stretch" style="width:100%;height:100%;" src="https://www.portlandmaps.com">
  <img src='https://jacksonvoelkel.com/geog575/images/backup_pdx_maps.png'>
</iframe>

## Open Source Stack

<img src='https://jacksonvoelkel.com/geog575/images/OS_GIS_Stack.jpg'>

##

<iframe class="stretch" style="width:100%;height:100%;" src="https://cmap.clackamas.us/maps/cmap">
  <img src='https://jacksonvoelkel.com/geog575/images/backup_c_maps.png'>
</iframe>

##

<br>

<h3><a href="https://latlon.io/apps/db_test/" target="_blank">Here's a simple example!</a></h3>

# Indexes

## Spatial Indexing

<br>

<img src='./img/bbox.png'>
<br><small><a href='https://postgis.net/workshops/postgis-intro/indexing.html' target='_blank'>"Spatial Indexing"</a></small>

## Indexes and PostGIS functions

<br>

Many of PostGIS's functions work naitively with spatial indexes, which add to it's incredible performance when compared to other tools.

# Big Take Aways

## Fast!

<br>

Not just in terms of time saved by being an all-in-one-shop, but the actual performance vs. proprietary systems!

## Simple!

<br>

No more moving data between applications, servers, and storage: just run analyses where the data lives!

## Great in Enterprise!

<br>

Allows teams to use the same data.

<br>

<p class='fragment'>This seems like such a 'duh' comment, but in the GIS world it's common to have multiple versions floating around multiple computers...</p>

## Accessible!

<br>

Though I didn't provide examples, in a large enterprise it is great to 'drive' PostGIS from anywhere!

<br>

<p class='fragment'>In fact, this entire presentation (including calling the queries and creating the maps) <a href='./postgis.Rmd' target="_blank">was made in R!</a></p>

## Expansive

<br>

- It is hard to give a talk on PostGIS and feel satisfied... I feel I just scratched the surface
- It does everything
- It connects to everything
- It's all in SQL
- It's shockingly fast
- It's free

##

# Thank You!

## Questions?

<br><br><br>

<a href="https://jacksonvoelkel.com" target="_blank">jacksonvoelkel.com</a>

<style type="text/css">code{white-space: pre-wrap;}</style>
<style type="text/css">.reveal section img {border: 0;}</style>
<script src="https://code.jquery.com/jquery-3.1.1.min.js" integrity="sha256-hVVnYaiADRTO2PzUGmuLJr8BLUSjGIZsDYGmIJLv2b8=" crossorigin="anonymous" onerror="rm()"></script>
<script src="https://ajax.googleapis.com/ajax/libs/jqueryui/1.10.4/jquery-ui.min.js" onerror="rm()"></script>
<link href="https://ajax.googleapis.com/ajax/libs/jqueryui/1.10.4/themes/smoothness/jquery-ui.min.css" rel="stylesheet" type="text/css" onerror="rm()"/>
<script>
$("p").has("img").each(function(){$(this).replaceWith($('<div class="stretch">' + this.innerHTML + '</div>'));});
$("img").each(function(){$(this).replaceWith($('<img src=' + this.src + '>' + this.innerHTML + '</img>'));});
$("li").addClass("fragment");
</script>