Skip to content

Instantly share code, notes, and snippets.

@worace

worace/gis_workshop.org

Last active February 6, 2023 02:56
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save worace/8c44fb2a514ddcb2ff98588813df7d02 to your computer and use it in GitHub Desktop.
Save worace/8c44fb2a514ddcb2ff98588813df7d02 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gis_workshop.org

GIS Intro

Spatial Stuff!

Existing Ecosystem:

“End-User” tools

  • Esri
  • ArcGIS
  • ???

GUIs, Shapefiles, etc

Database Integration

  • Indexing
  • Speed (GEOS-based)

PostGIS

Installing Postgres

Homebrew

brew install postgres
brew services start postgresql
createdb

#verify:
psql

Postgres.app

Not sure…download it from the website

Aside: Postgres Extensions

Installing Postgis (Postgres Extension)

Homebrew

brew install postgis

#test:
psql -c "create database postgis_test"
psql -d "postgis_test" -c "create extension postgis"
psql -c "drop database postgis_test"

Postgres.app

Should already be installed!

Postgis Basics

Data Types

  • Points
  • Polygons
  • LineStrings
  • Multi-things

Examples

-- Turing Point
select ST_Point(-104.996604, 39.750837);

-- Turing Polygon
select ST_GeometryType(ST_GeomFromText('POLYGON((-104.9964988231659 39.75115187915411,-104.99618232250214 39.75088379783092,-104.99657392501831 39.75055797451056,-104.99691188335419 39.7508260571017,-104.9964988231659 39.75115187915411))'));

-- How many points in the polygon?
select ST_NumPoints(ST_ExteriorRing(ST_GeomFromText('POLYGON((-104.9964988231659 39.75115187915411,-104.99618232250214 39.75088379783092,-104.99657392501831 39.75055797451056,-104.99691188335419 39.7508260571017,-104.9964988231659 39.75115187915411))')));

Aside: GIS Geometry Representations

What to polygons and linestrings look like?

Predicates / Query Functions

  • ST_Area
  • ST_Perimeter
  • ST_Distance
  • ST_Contains
  • ST_Intersects
  • ST_DWithin
select ST_Contains(ST_GeomFromText('POLYGON((-104.9964988231659 39.75115187915411,-104.99618232250214 39.75088379783092,-104.99657392501831 39.75055797451056,-104.99691188335419 39.7508260571017,-104.9964988231659 39.75115187915411))'),
                   ST_Point(-104.996604, 39.750837));

-- Union Station Point
select ST_Point(-105.000125, 39.753163);

select ST_Distance(ST_Point(-105.000125, 39.753163), ST_Point(-104.996604, 39.750837));
-- ASIDE: Units and Geography vs Geometry!?!?

select ST_Distance(ST_Point(-105.000125, 39.753163)::geography, ST_Point(-104.996604, 39.750837)::geography);

See the docs: https://postgis.net/docs/reference.html#Geometry_Constructors

Working with More Data

Data borrowed from this great tutorial: http://workshops.boundlessgeo.com/postgis-intro/index.html

Dir Setup and Data Bundle Download

# cd to wherever you put workshops and stuff
# cd ~/workshops
mkdir gis
cd gis

# download the data bundle
wget -O postgis_workshop.zip http://files.boundlessgeo.com/workshopmaterials/postgis-workshop-201401.zip

# unzip it
unzip postgis_workshop.zip

# take a look
tree postgis-workshop

# make a database for the workshop
psql -c "create database gis_workshop"
psql -d "gis_workshop" -c "create extension postgis"

# load some data

for f in ./postgis-workshop/data/*.shp
do
table_name=$(echo $f | sed 's/.*data\///' | sed 's/\.shp//')
echo "Loading $table_name from $f"
shp2pgsql -I -s 26918 $f $table_name | psql -d "gis_workshop" -U $(whoami)
done

psql -d "gis_workshop"

Examining Data

Tables:

gis_workshop=# \dt
               List of relations
 Schema |        Name         | Type  | Owner
--------+---------------------+-------+--------
 public | nyc_census_blocks   | table | horace
 public | nyc_homicides       | table | horace
 public | nyc_neighborhoods   | table | horace
 public | nyc_streets         | table | horace
 public | nyc_subway_stations | table | horace
 public | spatial_ref_sys     | table | horace

gis_workshop=# \d nyc_neighborhoods
                                     Table "public.nyc_neighborhoods"
  Column  |             Type             |                            Modifiers
----------+------------------------------+-----------------------------------------------------------------
 gid      | integer                      | not null default nextval('nyc_neighborhoods_gid_seq'::regclass)
 boroname | character varying(43)        |
 name     | character varying(64)        |
 geom     | geometry(MultiPolygon,26918) |
Indexes:
    "nyc_neighborhoods_pkey" PRIMARY KEY, btree (gid)
    "nyc_neighborhoods_geom_idx" gist (geom)

Crude Visualization

Aside: GIS Serialization Formats

  • WKT
  • WKB
  • GeoJSON

Aside: SRID’s

Consider:

select ST_AsGeoJSON(ST_SetSRID(geom, 4326)::geography) from nyc_neighborhoods limit 1;

Versus:

select ST_AsGeoJSON(ST_Transform(geom, 4326)::geography) from nyc_neighborhoods limit 1 offset 1;

Getting All of them:

psql -d gis_workshop -c "COPY (select ST_AsGeoJSON(ST_Transform(ST_Collect(ARRAY(select geom from nyc_neighborhoods limit 5)), 4326))) to STDOUT;" | pbcopy

Analysis

Aside: Power of GIS Joins

select count(*)
from nyc_homicides h
inner join nyc_neighborhoods n
on ST_Contains(n.geom, h.geom)
where n.boroname = 'Brooklyn';

What do you wanna know!?

Possible Q’s

  • Most populous borough?
  • Most populous borough by race?
  • # of homicides in most populous census district?
  • # of subway stations per borough? Per neighborhood?

GIS + Web

Take GIS Stuff from Database & Put it into a UI

GeoJSON

Maps

2 Problems

  • Getting Data to the map
  • Doing the Rendering / UI / etc
  • Write a lot of JS
  • Mapbox has SDKs for android and iOS as well!

Sinatra + Mapbox Skeleton

Gemfile

source 'https://rubygems.org'

gem 'sinatra'
gem 'sequel'
gem 'shotgun'
gem 'pg'

app.rb

require 'sinatra'

get '/' do
  File.read(File.join('public', 'index.html'))
end

public/index.html

<!DOCTYPE html>
<html lang="en">
  <head>
    <meta charset="utf-8">
    <title>GIS Workshop</title>
    <link rel="stylesheet" href="style.css">
    <script src="script.js"></script>
  </head>
  <body>
    <h1>Make a Map</h1>
    <div id='map' style="width: 800px; height: 800px;"></div>
  </body>
</html>

public/script.js

alert('hello');

Run It

shotgun app.rb

Mapbox Basics

Sign up and get a token

https://www.mapbox.com/studio/account/

Add JS Deps

No need for NPM, doing it old school.

<script src="https://code.jquery.com/jquery-3.2.1.min.js" integrity="sha256-hwg4gsxgFZhOsEEamdOYGBf13FyQuiTwlAQgxVSNgt4=" crossorigin="anonymous"></script>
<script src="https://api.mapbox.com/mapbox-gl-js/v0.39.1/mapbox-gl.js"></script>
<link href="https://api.mapbox.com/mapbox-gl-js/v0.39.1/mapbox-gl.css" rel='stylesheet' />

First Map

public/script.js:

const TOKEN = 'pk.eyJ1Ijoid29yYWNlIiwiYSI6ImNpeWMxOW1jcjAwYWUyd294ZzQ0YnMyZ3QifQ.ZaWekMcNTGFN-TmpPkf9AA';

mapboxgl.accessToken = TOKEN

function initMap() {
  return new mapboxgl.Map({
    container: 'map', // container id
    style: 'mapbox://styles/mapbox/streets-v9', // stylesheet location
    center: [-104.996604, 39.750837], // starting position [lng, lat]
    zoom: 13 // starting zoom
  });
}

$(document).ready(initMap);

Returning Neighborhood Data

app.rb

require 'sinatra'
require 'sequel'
require 'pg'
require 'json'

DB = Sequel.connect('postgres://@localhost:5432/gis_workshop')

get '/' do
  File.read(File.join('public', 'index.html'))
end

def neighborhoods_feature_collection
  query = 'select ST_AsGeoJSON(ST_Transform(geom, 4326)) as json from nyc_neighborhoods'
  rows = DB[query]
  features = rows.map do |r|
    JSON.parse(r[:json])
  end.map do |geom|
    {type: "Feature", geometry: geom}
  end
  {type: "FeatureCollection", features: features}
end

get '/neighborhoods' do
  content_type :json
  neighborhoods_feature_collection.to_json
end

script.js

const TOKEN = 'pk.eyJ1Ijoid29yYWNlIiwiYSI6ImNpeWMxOW1jcjAwYWUyd294ZzQ0YnMyZ3QifQ.ZaWekMcNTGFN-TmpPkf9AA';

mapboxgl.accessToken = TOKEN

console.log('running')

function initMap() {
  console.log('init');
  const map = new mapboxgl.Map({
    container: 'map', // container id
    style: 'mapbox://styles/mapbox/streets-v9', // stylesheet location
    center: [-73.991734, 40.698238], // starting position [lng, lat]
    zoom: 13 // starting zoom
  });

  map.on('load', () => {
    $.getJSON('/neighborhoods')
      .then(data => {
        console.log(data);
        map.addLayer({
          'id': 'neighborhoods',
          'type': 'fill',
          'source': {
            'type': 'geojson',
            'data': data
          },
          'layout': {},
          'paint': {
            'fill-color': '#088',
            'fill-opacity': 0.8
          }
        });

      });
  })
}

$(document).ready(initMap);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment