/*
Prereq:

Relational database experience, preferably PostgreSQL

SQL, SQL, SQL: If no have, focus on the functions

What is a spatial extension to a database?

What to take away: (Tools vs HD analogy)
1. Capability of PostGIS
2. Application to your own projects

Following along
*/

-- Create the Database

------ [Do following in another database.]
DROP DATABASE IF EXISTS workshop_1;
CREATE DATABASE workshop_1 ;

-- connect to new database in psql can use \connect workshop_1

-- Create Schemas [Explain use of schema]

CREATE SCHEMA postgis AUTHORIZATION postgres;
CREATE SCHEMA contrib AUTHORIZATION postgres;

-- Set Search Path [Explain]

ALTER DATABASE workshop_1 SET search_path = public, contrib, postgis;
CREATE EXTENSION postgis SCHEMA postgis;

------ [For variety, install hstore via PGAdmin otherwise run: CREATE EXTENSION hstore SCHEMA contrib;]

DROP TABLE IF EXISTS features;

CREATE TABLE features (id serial);
ALTER TABLE features ADD CONSTRAINT pk_features PRIMARY KEY (id);

-- PostGIS introduces new data types

ALTER TABLE features ADD COLUMN name varchar(20) NOT NULL;
ALTER TABLE features ADD COLUMN points geometry(point);
ALTER TABLE features ADD COLUMN linestrings geometry(linestring);
ALTER TABLE features ADD COLUMN polygons geometry(polygon);

CREATE INDEX ix_features_points ON features USING GiST (points);
CREATE INDEX ix_features_linestrings ON features USING GiST (linestrings);
CREATE INDEX ix_features_polygons ON features USING GiST (polygons);

------ [Modeling: Terrestrial, Political, Zoom levels]
------ [PostgreSQL does have native geometry types, but stay away from these.]

TRUNCATE TABLE features;

-- Points
INSERT INTO features (name, points) VALUES ('10 First Ave', ST_GeomFromText('POINT(4 10)'));
INSERT INTO features (name, points) VALUES ('15 First Ave', ST_GeomFromText('POINT(4 15)'));
INSERT INTO features (name, points) VALUES ('20 First Ave', ST_GeomFromText('POINT(4 20)'));

------ [Introduce OJ, not something we're formally teaching.]

-- OJ
SELECT name, points FROM features WHERE points IS NOT NULL
UNION ALL
SELECT name, linestrings FROM features WHERE linestrings IS NOT NULL
UNION ALL
SELECT name, polygons FROM features WHERE polygons IS NOT NULL;

-- Polygons [Draw polygon then translate into position.]
WITH 
    school As (SELECT 'POLYGON((0 0, 2 0, 2 1, 1 2, 0 1, 0 0))'::geometry(POLYGON) As geom),
	hospital As (SELECT 'POLYGON((0 1, 0 2, 1 2, 1 3, 2 3, 2 2, 3 2, 3 1, 2 1, 2 0, 1 0, 1 1, 0 1))'::geometry(POLYGON) As geom),
	town_hall As (SELECT 'POLYGON((0 0, 3 0, 3 3, 0 3, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1))'::geometry(POLYGON) as geom)
INSERT INTO features (name, polygons)
SELECT 'School' As name, ST_Translate(geom,22,24) FROM school
UNION ALL
SELECT 'Hospital' As name, ST_Translate(geom,0,0) FROM hospital
UNION ALL
SELECT 'Town Hall' As name, ST_Translate(geom,14,14) FROM town_hall;

------ [Note parenthesis]

-- Linestrings

INSERT INTO features (name, linestrings) VALUES ('Main St', ST_GeomFromText('LINESTRING(0 4, 14 4, 14 24, 28 24)'));
INSERT INTO features (name, linestrings) VALUES ('First Ave', ST_GeomFromText('LINESTRING(4 0, 4 28)'));
INSERT INTO features (name, linestrings) VALUES ('Second Ave', ST_GeomFromText('LINESTRING(8 0, 8 28)'));

INSERT INTO features (name, linestrings) 
VALUES (
    'ER', 
    ST_MakeLine(
        (SELECT ST_ClosestPoint((SELECT linestrings FROM features WHERE name = 'Main St'),(SELECT polygons FROM features WHERE name = 'Hospital'))),
        (SELECT ST_ClosestPoint((SELECT polygons FROM features WHERE name = 'Hospital'),(SELECT linestrings FROM features WHERE name = 'Main St')))
    )
);
        
------ [Don't use ST_PointFromText, ST_LineStringFromText, ST_PolygonFromText]
------ [Use always ST_GeomFromText, ST_Point or ST_MakePoint]

------ [Show raw table data. Note that you can't make sense out of the geometry columns.]
------ [Alternative table designs. All sub-types in one column. All in separate tables.]

------ [Functions: Need to know only a few functions, use intuition on which function works on which geometry.]
------ [Some functions are heavily overloaded. ST_AsGeoJSON, one that takes gometry, one that takes geography, one that specifies significant digits.]
------ [Different versions of same function exists depending on how PostGIS was compiled.]
------ [Function always prefixed with ST, to distinguish from other PostgreSQL functions, internal functions.]

SELECT name, ST_Length(linestrings) FROM features WHERE linestrings IS NOT NULL;
SELECT name, ST_Length(points) FROM features WHERE points IS NOT NULL;
SELECT name, ST_Length(polygons) FROM features WHERE polygons IS NOT NULL;

SELECT name, ST_Area(polygons) FROM features WHERE polygons IS NOT NULL;
SELECT name, ST_Perimeter(polygons) FROM features WHERE polygons IS NOT NULL;

SELECT name, ST_NPoints(polygons) FROM features WHERE polygons IS NOT NULL;
SELECT name, ST_NPoints(linestrings) FROM features WHERE linestrings IS NOT NULL; -- [SELECT name, ST_NumPoints(linestrings) FROM features WHERE linestrings IS NOT NULL;]

SELECT name, ST_AsText(ST_Centroid(polygons)) FROM features WHERE polygons IS NOT NULL;
SELECT name, ST_AsText(ST_Centroid(linestrings)) FROM features WHERE linestrings IS NOT NULL;

SELECT name, ST_XMax(linestrings) FROM features WHERE linestrings IS NOT NULL;
SELECT name, ST_XMin(linestrings) FROM features WHERE linestrings IS NOT NULL;
SELECT ST_X(linestrings) FROM features; -- [This doesn't make sense.]

-- Output Functions

SELECT id, name, ST_AsText(points), ST_AsText(linestrings), ST_AsText(polygons) FROM features; -- [ST_AsWKB(), ST_AsKML, ST_GML, ST_AsGeoJSON are some other popular output functions.]

-- ST_Intersects

-- Point and LineString
WITH 
    A AS (SELECT name, points As geom FROM features WHERE points IS NOT NULL),
    B AS (SELECT name, linestrings As geom FROM features WHERE linestrings IS NOT NULL)
SELECT A.name, B.name, ST_Intersects(A.geom,B.geom) FROM A CROSS JOIN B;

-- Point and Polygon
WITH 
    A AS (SELECT name, points As geom FROM features WHERE points IS NOT NULL),
    B AS (SELECT name, polygons As geom FROM features WHERE polygons IS NOT NULL)
SELECT A.name, B.name, ST_Intersects(A.geom,B.geom), ST_Within(A.geom,B.geom) FROM A CROSS JOIN B;

-- LineString and Polygon
WITH 
    A AS (SELECT name, linestrings As geom FROM features WHERE linestrings IS NOT NULL),
    B AS (SELECT name, polygons As geom FROM features WHERE polygons IS NOT NULL)
SELECT A.name, B.name, ST_Intersects(A.geom,B.geom) FROM A CROSS JOIN B;

-- LineString and LineString
WITH 
    A AS (SELECT name, linestrings As geom FROM features WHERE linestrings IS NOT NULL),
    B AS (SELECT name, linestrings As geom FROM features WHERE linestrings IS NOT NULL)
SELECT A.name, B.name, ST_Intersects(A.geom,B.geom), ST_Touches(A.geom,B.geom), ST_Crosses(A.geom,B.geom) FROM A CROSS JOIN B;

DROP TABLE IF EXISTS regions;
CREATE TABLE regions (id serial);
ALTER TABLE regions ADD CONSTRAINT pk_regions PRIMARY KEY (id);
ALTER TABLE regions ADD COLUMN name varchar(20) NOT NULL;
ALTER TABLE regions ADD COLUMN polygons geometry(polygon);
CREATE INDEX ix_regions_polygons ON regions USING GiST (polygons);


INSERT INTO regions (name, polygons) SELECT 'NW' As name, ST_GeomFromText('POLYGON((0 14, 14 14, 14 28, 0 28, 0 14))');
INSERT INTO regions (name, polygons) SELECT 'NE' As name, ST_GeomFromText('POLYGON((14 14, 28 14, 28 28, 14 28, 14 14))');
INSERT INTO regions (name, polygons) SELECT 'SE' As name, ST_GeomFromText('POLYGON((14 0, 28 0, 28 14, 14 14, 14 0))');
INSERT INTO regions (name, polygons) SELECT 'SW' As name, ST_GeomFromText('POLYGON((0 0, 14 0, 14 14, 0 14, 0 0))');
INSERT INTO regions (name, polygons) SELECT 'Town' As name, ST_GeomFromText('POLYGON((0 0, 28 0, 28 28, 0 28, 0 0))');

-- Polygon and polygon

WITH 
    A AS (SELECT name, polygons As geom FROM regions),
    B AS (SELECT name, polygons As geom FROM features WHERE polygons IS NOT NULL)
SELECT A.name, B.name, ST_Intersects(A.geom,B.geom), ST_Overlaps(A.geom,B.geom), ST_Contains(A.geom,B.geom), ST_Contains(B.geom,A.geom) FROM A CROSS JOIN B
ORDER BY B.name, A.name;

-- Intersections

WITH 
    A AS (SELECT name, linestrings As geom FROM features WHERE name = 'First Ave'),
    B AS (SELECT name, linestrings As geom FROM features WHERE name = 'Main St')
SELECT A.name, B.name, ST_AsText(ST_Intersection(A.geom,B.geom)) FROM A CROSS JOIN B;
 
WITH 
    A AS (SELECT name, polygons As geom FROM regions WHERE name = 'NE'),
    B AS (SELECT name, polygons As geom FROM regions WHERE name = 'Town')
SELECT A.name, B.name, ST_AsText(ST_Intersection(A.geom,B.geom)) FROM A CROSS JOIN B;

WITH 
    A AS (SELECT name, polygons As geom FROM regions WHERE name = 'NE'),
    B AS (SELECT name, polygons As geom FROM regions WHERE name = 'SW')
SELECT A.name, B.name, ST_AsText(ST_Intersection(A.geom,B.geom)) FROM A CROSS JOIN B;

WITH 
    A AS (SELECT name, polygons As geom FROM regions WHERE name = 'NE'),
    B AS (SELECT name, polygons As geom FROM regions WHERE name = 'NW')
SELECT A.name, B.name, ST_AsText(ST_Intersection(A.geom,B.geom)) FROM A CROSS JOIN B;

WITH 
    A AS (SELECT name, polygons As geom FROM regions WHERE name = 'NE'),
    B AS (SELECT name, polygons As geom FROM features WHERE name = 'Hospital')
SELECT A.name, B.name, ST_AsText(ST_Intersection(A.geom,B.geom)) FROM A CROSS JOIN B;

-- Buffer

DROP TABLE IF EXISTS buffers;

CREATE TABLE buffers (
    id serial NOT NULL, 
    name varchar(20) NOT NULL, 
    polygons geometry(polygon) NOT NULL, 
    CONSTRAINT pk_buffers PRIMARY KEY (id)
);

CREATE INDEX ix_buffers_polygons ON regions USING GiST (polygons);

INSERT INTO buffers (name, polygons) 
SELECT 'School Buffer Zone', ST_Buffer(polygons,.5) FROM features WHERE name = 'School'
UNION ALL
SELECT 'Hospital Buffer Zone', ST_Buffer(polygons,.5) FROM features WHERE name = 'Hospital'
UNION ALL
SELECT 'Main St Buffer Zone', ST_Buffer(linestrings,.5) FROM features WHERE name = 'Main St';

-- Equality

INSERT INTO features (name, linestrings) VALUES ('E Solidus Ln', ST_GeomFromText('LINESTRING(4 10, 8 14)'));
INSERT INTO features (name, linestrings) VALUES ('W Solidus Ln', ST_GeomFromText('LINESTRING(4 14, 8 10)'));

WITH 
    A AS (SELECT name, linestrings As geom FROM features WHERE name = 'E Solidus Ln'),
    B AS (SELECT name, linestrings As geom FROM features WHERE name = 'W Solidus Ln')
SELECT A.geom = B.geom FROM A CROSS JOIN B;

WITH 
    A AS (SELECT name, linestrings As geom FROM features WHERE name = 'E Solidus Ln'),
    B AS (SELECT name, linestrings As geom FROM features WHERE name = 'W Solidus Ln')
SELECT ST_Equals(A.geom,B.geom) FROM A CROSS JOIN B;

SELECT Box2D(linestrings) FROM features WHERE name LIKE '%Solidus Ln'; -- [BOX is not a geometry!]

-- Distance
------ [Distance is always the closest points between two geometries.]

WITH X As (SELECT name, points AS geom FROM features WHERE points IS NOT NULL)
SELECT L.name, R.name, ST_Distance(L.geom,R.geom) 
FROM X L INNER JOIN X R ON L.name <> R.name;

WITH X As (SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL)
SELECT L.name, R.name, ST_Distance(L.geom,R.geom) 
FROM X L INNER JOIN X R ON NOT ST_Equals(L.geom,R.geom);

WITH X As (SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL)
SELECT L.name, R.name, ST_Distance(L.geom,R.geom) 
FROM X L INNER JOIN X R ON NOT ST_Equals(L.geom,R.geom);

SELECT ST_AsText(ST_Centroid(polygons)) FROM regions WHERE name = 'Town'; -- [Centroid of the Town]

WITH 
    X As (
        SELECT name, points AS geom FROM features WHERE points IS NOT NULL UNION ALL 
        SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL UNION ALL 
        SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
    )
SELECT L.name, R.name, ST_Distance(L.geom,R.geom) 
FROM X L INNER JOIN X R ON NOT ST_Equals(L.geom,R.geom);

-- Within X distance

WITH 
    X As (
        SELECT name, points AS geom FROM features WHERE points IS NOT NULL UNION ALL 
        SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL UNION ALL 
        SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
    ),
    Y AS (SELECT ST_Centroid(polygons) AS geom FROM regions WHERE name = 'Town')
SELECT name FROM X, Y WHERE ST_Distance(X.geom,Y.geom) <= 7.5;

WITH 
    X As (
        SELECT name, points AS geom FROM features WHERE points IS NOT NULL UNION ALL 
        SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL UNION ALL 
        SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
    ),
    Y AS (SELECT ST_Centroid(polygons) AS geom FROM regions WHERE name = 'Town')
SELECT name FROM X, Y WHERE ST_Intersects(ST_Buffer(Y.geom,7.5),X.geom);

WITH 
    X As (
        SELECT name, points AS geom FROM features WHERE points IS NOT NULL UNION ALL 
        SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL UNION ALL 
        SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
    ),
    Y AS (SELECT ST_Centroid(polygons) AS geom FROM regions WHERE name = 'Town')
SELECT name FROM X, Y WHERE ST_DWithin(Y.geom,X.geom,7.5);

---- KNN Bounding Box

WITH 
    X As (
        SELECT name, points AS geom FROM features WHERE points IS NOT NULL UNION ALL 
        SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL UNION ALL 
        SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
    ),
    Y AS (SELECT ST_Centroid(polygons) AS geom FROM regions WHERE name = 'Town')
SELECT name 
FROM X, Y
ORDER BY X.geom <#> Y.geom;
------ [This will NOT use a spatial index. Not good.]

WITH 
    X As (
        SELECT name, points AS geom FROM features WHERE points IS NOT NULL UNION ALL 
        SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL UNION ALL 
        SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
    )
SELECT name 
FROM X
ORDER BY X.geom <#> ST_GeomFromText('POINT(14 14)');
------ [This will use a spatial index. Better.]

---- KNN Bounding Box
WITH 
    X As (
        SELECT name, points AS geom FROM features WHERE points IS NOT NULL UNION ALL 
        SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL UNION ALL 
        SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
    )
SELECT name 
FROM X
ORDER BY X.geom <#> ST_GeomFromText('POINT(1 15)');

---- KNN Centroid
WITH 
    X As (
        SELECT name, points AS geom FROM features WHERE points IS NOT NULL UNION ALL 
        SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL UNION ALL 
        SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
    )
SELECT name 
FROM X
ORDER BY X.geom <-> ST_GeomFromText('POINT(1 15)');

-- MULTI

SELECT ST_AsText(ST_Collect(points)) FROM features;
SELECT ST_AsText(ST_Collect(linstrsings)) FROM features;
SELECT ST_AsText(ST_Collect(polygons)) FROM features;
SELECT ST_AsText(polygons) FROM features WHERE polygons IS NOT NULL;
SELECT ST_AsText(ST_Multi(polygons)) FROM features WHERE polygons IS NOT NULL;

WITH 
    X As (
        SELECT name, points AS geom FROM features WHERE points IS NOT NULL UNION ALL 
        SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL UNION ALL 
        SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
    )
SELECT ST_AsText(ST_Collect(geom)) FROM X;

SELECT name, ST_NumGeometries(geom), geom FROM data.us ORDER BY name; -- [US Bureau of Census]

------ [SRID encapsulates Geoid, Datum, Projection]

SELECT srid, srtext::varchar(50) 
FROM spatial_ref_sys 
WHERE srtext ILIKE '%83%' AND srtext ILIKE '%California%' AND proj4text ILIKE '%ft%';

DROP TABLE IF EXISTS ni; 

SELECT name, geom INTO ni 
FROM (
    SELECT name, points AS geom FROM features WHERE points IS NOT NULL
    UNION ALL
    SELECT name, linestrings AS geom FROM features WHERE linestrings IS NOT NULL
    UNION ALL
    SELECT name, polygons AS geom FROM features WHERE polygons IS NOT NULL
) X

ALTER TABLE ni ALTER COLUMN geom TYPE geometry(geometry,4326) USING ST_SetSRID(geom,4326);
CREATE INDEX ix_ni_geom ON ni USING GiST (geom);

------ [ST_SetSRID only resets the meta data. Apply ST_Transform to do the actual transformation]
------ [Can't use ST_Transform from NO SRID to something with SRID.]

SELECT ST_Transform(points,4326) FROM features WHERE points IS NOT NULL;
SELECT ST_SetSRID(points,9999) FROM features WHERE points IS NOT NULL;
SELECT ST_Transform(ST_SetSRID(points,3857),4326) FROM features WHERE points IS NOT NULL; -- [Mercator, used by GoogleMaps]

-- NULL Island [http://www.nullisland.com/geography.html]

SELECT name, ST_Intersection(geom,ST_SetSRID('BOX(0 0, 28 28)'::Box2D,4326)) 
FROM data.countries 
WHERE geom && 'BOX(0 0, 28 28)'::Box2D; -- [OJ]

SELECT geom FROM ni; -- [OJ]

-- Geography

ALTER TABLE ni ADD COLUMN geog geography;

UPDATE ni SET geog = geom::geography; -- [SRID has to be in a Lon-Lat. This won't work: UPDATE ni SET geog = ST_Transform(geom,3857)::geography;]

ALTER TABLE data.countries DROP COLUMN IF EXISTS geog;
ALTER TABLE data.countries ADD COLUMN geog geography;
UPDATE data.countries SET geog = geom::geography;

WITH 
    A AS (SELECT geom FROM data.countries WHERE name = 'Japan'),
    B AS (SELECT geom FROM ni WHERE name = 'Town Hall')
SELECT ROUND(ST_Distance(A.geom,B.geom)) FROM A CROSS JOIN B;

WITH 
    A AS (SELECT geom FROM data.countries WHERE name = 'Japan'),
    B AS (SELECT geom FROM ni WHERE name = 'Town Hall')
SELECT ROUND(ST_Distance_Sphere(A.geom,B.geom) * .000621371) AS miles FROM A CROSS JOIN B; -- [Assumes earth to be a perfect sphere.]

WITH 
    A AS (SELECT geom FROM data.countries WHERE name = 'Japan'),
    B AS (SELECT geom FROM ni WHERE name = 'Town Hall')
SELECT ROUND(ST_Distance_Spheroid(A.geom,B.geom,'SPHEROID["WGS 84",6378137,298.257223563]'::spheroid) * .000621371) AS miles FROM A CROSS JOIN B; 

WITH 
    A AS (SELECT geog FROM data.countries WHERE name = 'Japan'),
    B AS (SELECT geog FROM ni WHERE name = 'Town Hall')
SELECT ROUND(ST_Distance(A.geog,B.geog) * .000621371) AS miles FROM A CROSS JOIN B;

------ [ITA: Import cities with GUI, rename imported to world_cities, move to data schema using ALTER world_cities SET SCHEMA data;]

WITH 
    A AS (SELECT name, geog::geometry(POINT,4326) As geom FROM data.world_cities WHERE name = 'Houston'),
    B As (SELECT name, geog::geometry(POINT,4326) As geom FROM data.world_cities WHERE name = 'Mumbai')
SELECT ST_Segmentize(ST_MakeLine(A.geom,B.geom)::geography, 100)::geometry FROM A CROSS JOIN B;

WITH 
    A AS (SELECT name, geog::geometry(POINT,4326) As geom FROM data.world_cities WHERE name = 'Houston'),
    B As (SELECT name, geog::geometry(POINT,4326) As geom FROM data.world_cities WHERE name = 'Mumbai')
SELECT ST_MakeLine(A.geom,B.geom) FROM A CROSS JOIN B;

SELECT ROUND(ST_Area(geog)/1000000) FROM data.countries WHERE name = 'United States';
SELECT ROUND(ST_Area(geom)) FROM data.countries WHERE name = 'United States';

/* Importing BART data using shp2pgsql
"C:\Program Files\PostgreSQL\9.4\bin\shp2pgsql" -D -I -t 2D -s 4269 C:\Workshop_1\Data\Bart\BART_13 data.bart_lines | "C:\Program Files\PostgreSQL\9.4\bin\psql" -p 5433 -U postgres -d workshop_1
"C:\Program Files\PostgreSQL\9.4\bin\shp2pgsql" -D -I -t 2D -s 4269 C:\Workshop_1\Data\Bart\BART_Sta_13 data.bart_stations | "C:\Program Files\PostgreSQL\9.4\bin\psql" -p 5433 -U postgres -d workshop_1
*/

------ [Data is in SRID 4269 (NAD83 Lon Lat)]
SELECT * FROM data.bart_lines;
SELECT line, ST_Length(geom) FROM data.bart_lines; -- [This would give degree measure.]
SELECT line, ST_Length(ST_Transform(geom,2163)) * .000621371 FROM data.bart_lines; -- [Transform to National Atlas Equatorial Something (meter-based), equal-area projection, distance perserving.]

-- Station within 10 KM [2163 National Atlas (Distance Preserving)]
WITH X As (SELECT ST_Transform(ST_GeomFromText('POINT(-122.365376 37.593817)',4269),2163) As geom)
SELECT station, ST_Distance(ST_Transform(S.geom,2163), X.geom)::integer/1000 As dist_mtr
FROM data.bart_stations S INNER JOIN X
ON ST_DWithin(ST_Transform(S.geom,2163), X.geom, 10000)

-- Three closest train stops using centroid distance
WITH X As (SELECT ST_GeomFromText('POINT(-122.365376 37.593817)',4269) As geom)
SELECT 
    station, 
    geom <-> (SELECT geom FROM X) As dist_knn 
FROM data.bart_stations
ORDER BY geom <-> (SELECT geom FROM X) LIMIT 3;

-- Geocoding

---- Tiger normalization

------ [ftp://ftp2.census.gov/geo/tiger/TIGER2014/STATE]
CREATE EXTENSION fuzzystrmatch SCHEMA contrib;

--disconnect database and reconnect (so get new alter database search_path)
CREATE EXTENSION postgis_tiger_geocoder; -- [This creates two new schemas automatically.]
CREATE EXTENSION address_standardizer SCHEMA contrib;

SELECT 'Non-PAGC', address, streetname, streettypeabbrev, zip FROM normalize_address('ONE COURT Place, Boston, MA 02109')
UNION ALL
SELECT 'PAGC', address, streetname, streettypeabbrev, zip FROM pagc_normalize_address('ONE COURT Place, Boston, MA 02109');


-- for this part need to load already prepped tiger_data_chic.sql.bz2 (refer to load script)
---- Geocoding

SELECT * FROM geocode('333 North Dearborn Street, Chicago, IL 60654',1);
SELECT * FROM geocode('6666 North Dearborn Street, Chicago, IL 60654',1);
SELECT * FROM geocode('999 Darth Vader Street, Chicago, IL 60654',1);
SELECT * FROM geocode('10 North, Chicago, IL 60654',10);

SELECT pprint_addy(addy), ST_AsText(geomout) FROM geocode('333 North Dearborn Street, Chicago, IL 60654',1);
SELECT pprint_addy(addy), ST_AsText(ST_SnapToGrid(geomout,0.001)) FROM geocode('333 North Dearborn Street, Chicago, IL 60654',1);
SELECT pprint_addy(addy), ST_AsText(ST_SnapToGrid(geomout,0.001)), rating FROM geocode_intersection('Lake', 'Clark', 'IL', 'Chicago');

------ [Reverse geocoding is also available, but outside the scope of this workshop.]

-- PGRouting
------ [Map of London Tube: http://www.tfl.gov.uk/maps/track]
------ [Tube lines imported using shp2pgsql. Stations is a CSV file, imported using psql.]
------ [Cleanup raw data which contains superfluous Z and M coordinates: ALTER TABLE data.london_tube_lines ALTER COLUMN geom TYPE geometry(LineString, 4326) USING ST_Force2D(geom);]
CREATE EXTENSION pgRouting;
ALTER EXTENSION pgRouting SET SCHEMA contrib;

CREATE TABLE data.london_stations (
   station character varying(150), 
   os_x integer, 
   os_y integer, 
   latitude numeric(10,6),
   longitude numeric(10,6),
   zone varchar(20), 
   postcode varchar(50),
   CONSTRAINT pk_london_stations PRIMARY KEY (station) 
);

---- Load data using psql 
\cd C:/Workshop_1/Data/London; -- [Note the Unix-style slashes]
\copy data.london_stations FROM 'London stations.csv' DELIMITER ',' CSV HEADER;

ALTER TABLE data.london_stations DROP COLUMN IF EXISTS geom;
ALTER TABLE data.london_stations ADD COLUMN geom geometry(POINT,4326);

DELETE FROM data.london_stations WHERE station = 'Euston Square';
UPDATE data.london_stations SET geom = ST_SetSRID(ST_Point(longitude,latitude),4326);

------ [Stations imported simply to name stations. The verticies are already in the london_tube_lines table.]
ALTER TABLE data.london_tube_lines DROP COLUMN IF EXISTS "source";
ALTER TABLE data.london_tube_lines DROP COLUMN IF EXISTS "target";

ALTER TABLE data.london_tube_lines ADD COLUMN "source" integer;
ALTER TABLE data.london_tube_lines ADD COLUMN "target" integer;

------ [The next super function populates source and target, creates london_tube_lines_vertices_pgr table.]
SELECT pgr_createTopology('data.london_tube_lines', 0.00001, 'geom', 'gid');

SELECT * FROM data.london_tube_lines_vertices_pgr;

ALTER TABLE data.london_stations DROP COLUMN IF EXISTS station_id;
ALTER TABLE data.london_stations ADD COLUMN station_id integer;

UPDATE data.london_stations 
SET station_id = X.id
FROM data.london_tube_lines_vertices_pgr X
WHERE ST_DWithin(data.london_stations.geom, X.the_geom, 0.000001);

DELETE FROM data.london_stations WHERE station_id IS NULL;
------ [ITA: SELECT count(*) As cnt, count(station_id) As tot_assigned, count(distinct station_id) As tot_distinct FROM data.london_stations;]

CREATE INDEX ix_source_london_tube_lines ON data.london_tube_lines(source);
CREATE INDEX ix_target_london_tube_lines ON data.london_tube_lines(target);

ALTER TABLE data.london_tube_lines DROP COLUMN IF EXISTS length;
ALTER TABLE data.london_tube_lines ADD COLUMN length float8;
UPDATE data.london_tube_lines SET length = ST_Length(geom::geography);

SELECT gid, name, length, source, Y.station, target, Z.station  
FROM 
    data.london_tube_lines X INNER JOIN 
    data.london_stations Y ON X.source = Y.station_id INNER JOIN 
    data.london_stations Z ON X.target = Z.station_id;

SELECT seq, id1 AS node, id2 AS edge, cost 
FROM 
    pgr_dijkstra('
        SELECT gid AS id, source::int4, target::int4, length::float8 AS cost FROM data.london_tube_lines',
        (SELECT station_id FROM data.london_stations WHERE station = 'Chesham'),
        (SELECT station_id FROM data.london_stations WHERE station = 'West Croydon'), 
        false, 
        false
    );

SELECT geom FROM data.london_tube_lines
WHERE gid in (
    SELECT id2
        FROM 
            pgr_dijkstra('
                SELECT gid AS id, source::int4, target::int4, length::float8 AS cost FROM data.london_tube_lines',
                (SELECT station_id FROM data.london_stations WHERE station = 'Chesham'),
                (SELECT station_id FROM data.london_stations WHERE station = 'West Croydon'), 
            false, 
            false
    )
);


SELECT seq, S.station, L.name, cost * .000621371 AS miles
FROM 
    pgr_dijkstra('
        SELECT gid AS id, source::int4, target::int4, length::float8 AS cost FROM data.london_tube_lines',
        (SELECT station_id FROM data.london_stations WHERE station = 'Chesham'),
        (SELECT station_id FROM data.london_stations WHERE station = 'West Croydon'), 
        false, 
        false
    ) R
    INNER JOIN data.london_stations S ON R.id1 = S.station_id
    LEFT JOIN data.london_tube_lines L ON R.id2 = L.gid
ORDER BY R.seq;

-- 3D

CREATE EXTENSION postgis_sfcgal SCHEMA postgis; -- this is only in PostGIS 2.2.  For 2.1m need to run the sfcgal.sql script 

--SFCGAL is not needed for these examples -- only needed for ST_Extrude we will show later
---- 3D Linestring
SELECT '<shape><appearance><material ambientintensity="0.200" containerfield="material" shininess="0.25" emissivecolor="1 .25 .25"></material></appearance>' || ST_AsX3D('LINESTRING Z(100 200 0, 200 200 0, 200 300 0, 200 350 0, 200 350 10, 200 350 50)'::geometry) ||  '</shape>' As x3d;
                       
---- Closed PolyhedralSurface with X3D
SELECT '<Shape><Appearance><Material ambientIntensity="0.200" containerField="material" shininess="0.200" diffuseColor="1 0 0"/></Appearance>' || 
    ST_AsX3D(
        'POLYHEDRALSURFACE Z( 
            ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)), 
            ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)), 
            ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)), 
            ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)), 
            ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0)), 
            ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1))
        )'::geometry
    ) || 
    '</Shape>';
                    
---- Closed TIN with X3D
SELECT 
    '<Shape><Appearance><Material ambientIntensity="0.200" containerField="material" shininess="0.200" diffuseColor="0.8 0 0.6"/></Appearance>' || 
    replace(
        ST_AsX3D(
            'TIN(
                ((0 0 0, 1 0 0, 0 1 0, 0 0 0)),
                ((0 0 0, 1 0 0, 0 0 1, 0 0 0)),
                ((0 0 0, 0 0 1, 0 1 0, 0 0 0)),
                ((0 0 1, 1 0 0, 0 1 0, 0 0 1))
            )'::geometry),  
        '&lt;IndexedTriangleSet',
        '&lt;IndexedTriangleSet solid="false" ') || 
    '</Shape>';     