-- Installing PostGIS
CREATE DATABASE gisdb;
\c gisdb -- if in psql

ALTER DATABASE gisdb SET search_path=public,postgis,tiger,contrib;
\c gisdb -- if in psql
CREATE SCHEMA IF NOT EXISTS postgis;
CREATE SCHEMA IF NOT EXISTS contrib;
CREATE EXTENSION postgis SCHEMA postgis; --- geometry, geography data types
-- separated in PostGIS 3.0 as separate extension
CREATE EXTENSION postgis_raster SCHEMA postgis; -- used to be part of postgis extension
CREATE EXTENSION postgis_sfcgal SCHEMA postgis; -- much of the functionality like dealing with 2d TINS and Polyhedralsurfaces now supported in PostGIS proper
CREATE EXTENSION postgis_topology;
-- used by postgis_tiger_geocoder
CREATE EXTENSION fuzzystrmatch SCHEMA contrib; -- note part of PostGIS
CREATE EXTENSION address_standardizer SCHEMA contrib; -- for standardizing addresses
CREATE EXTENSION postgis_tiger_geocoder; -- needs postgis and fuzzystrmatch, can use address_standardizer

-- How you upgrade 
SELECT postgis_extensions_upgrade();

-- Verify version of PostGIS
SELECT postgis_full_version();

/**POSTGIS="3.2.0dev 3.1.0rc1-207-g296e1bfdb" [EXTENSION] PGSQL="140" GEOS="3.10.0dev-CAPI-1.15.0" SFCGAL="1.3.8" PROJ="7.1.1" GDAL="GDAL 3.2.3, released 2021/04/27" LIBXML="2.9.9" LIBJSON="0.12" LIBPROTOBUF="1.2.1" WAGYU="0.5.0 (Internal)" TOPOLOGY RASTER**/

-- useful for linking to foreign data
CREATE EXTENSION ogr_fdw SCHEMA postgis;

-- INSTALL http extension (not part of postgis) - https://github.com/pramsey/pgsql-http 
CREATE EXTENSION http SCHEMA contrib;

SELECT http_set_curlopt('CURLOPT_SSL_VERIFYPEER','0');

-- SLIDE: extract features from a geojson collection
DROP TABLE IF EXISTS transportation_districts;
WITH d AS (SELECT  h.content::jsonb AS content
FROM http_get('https://bostonopendata-boston.opendata.arcgis.com/datasets/f953f854b672496fb0aa18ad92278f07_5.geojson?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D') AS h)
SELECT
    r->'properties'->>'OBJECTID_1' AS object_id,
    r->'properties'->>'DIST' AS district, 
    ST_GeomFromGeoJSON(r->>'geometry')::geography AS geog,
    r AS data_raw
INTO transportation_districts
FROM d, jsonb_array_elements(d.content->'features') AS r;

-- SLIDE: Visualize and measure
SELECT district, geog, ST_Area(geog) AS sq_meters
FROM transportation_districts;

-- SLIDE: Create points from long /lat - 
-- 19477 rows
-- building permits - 17631 rows
DROP TABLE IF EXISTS recent_building_permits;
SELECT ST_Point((r->>'long')::float8, (r->>'lat')::float8)::geography(POINT, 4326) AS geog,
	r->>'permitnumber' AS permit_number, (r->>'issued_date')::date AS issued_date, 
	r->>'worktype' AS work_type, r->>'description' AS description,
	(r->>'total_fees')::numeric AS total_fees
INTO recent_building_permits
FROM http_get('https://data.boston.gov/api/3/action/datastore_search_sql?sql=' 
	|| urlencode($sql$SELECT lat,long, permitnumber,issued_date,worktype,description,total_fees 
	FROM "6ddcd912-32a0-43df-9908-63574f8c7e77" WHERE "issued_date" > '2021-01-01' $sql$) ) AS h,
	jsonb_array_elements(h.content::jsonb->'result'->'records') AS r;

-- SLIDE: Visualize building permits
SELECT geog, work_type FROM recent_building_permits;

-- SLIDE: indexes make things faster
CREATE INDEX ix_recent_building_permits_geog 
    ON recent_building_permits USING gist(geog);
    
CREATE INDEX ix_transportation_districts_geog 
    ON  transportation_districts USING gist(geog);
 
CREATE INDEX ix_recent_building_permits_work_type
    ON recent_building_permits USING btree(work_type);

-- SLIDE: How much plumbing is going on in each district
SELECT d.district, d.geog, COUNT(p.work_type) 
FROM transportation_districts AS d 
    INNER JOIN recent_building_permits AS p ON ST_Covers(d.geog, p.geog)
WHERE work_type = 'PLUMBING'
GROUP BY d.district, d.geog
ORDER BY d.district, d.geog;

-- SLIDE : Geocoding with PostGIS tiger geocoder 
-- detailed load process http://postgis.net/docs/postgis_installation.html#tiger_geocoder_loading_data
SELECT install_missing_indexes();
CREATE EXTENSION postgis_tiger_geocoder;
-- SLIDE: Basic geocode
-- tiger data is in north american datum which is pretty close to the standard WGS 84 long lat
-- pgAdmin doesn't overlay a map if data is not in WGS 84 (4326)
SELECT ST_Buffer(ST_SetSRID(geomout, 4326)::geography,2), (addy).*
FROM geocode('1 Devonshire Place, PH 301, Boston, MA 02109',1);

SELECT ST_AsText(ST_SetSRID(geomout, 4326)) , (addy).*
FROM geocode('1 Devonshire Place, PH 301, Boston, MA 02109',1);

SELECT ST_Buffer(ST_SetSRID(geomout, 4326),0.00003,1), (addy).*
FROM geocode('75 State St, Boston, MA',1);

-- SLIDE: Geocoding table of data
-- women owned businesses in Boston Area - https://data.boston.gov/dataset/women-owned-businesses/resource/91ff216b-4f54-4a2c-af00-2c0f987070ee

-- 13 rows
DROP TABLE IF EXISTS women_owned;
SELECT (r->>'_id')::bigint AS id, r::jsonb AS data_raw
INTO women_owned
FROM http_get('https://data.boston.gov/api/3/action/datastore_search_sql?sql=' 
	|| urlencode($sql$SELECT * 
	FROM "91ff216b-4f54-4a2c-af00-2c0f987070ee" WHERE "Business Type" = 'Education' $sql$) ) AS h,
	jsonb_array_elements(h.content::jsonb->'result'->'records') AS r;

SELECT data_raw FROM women_owned;
-- 824 ms
DROP TABLE IF EXISTS women_owned_geocoded;
WITH wo AS (SELECT id, data_raw->>'Business Name' As business_name,
    data_raw->>'Physical Location/Address' AS address_1, 
    data_raw->>'Business Zipcode' AS zip_code, data_raw
            FROM women_owned)
SELECT wo.id, wo.business_name, 
    wo.address_1 || ',Boston, MA ' || wo.zip_code AS orig_full_address,
    ST_SetSRID(g.geomout,4326)::geography AS geog, 
    (g.addy).*, g.rating, wo.data_raw
    INTO women_owned_geocoded
FROM  wo LEFT JOIN LATERAL geocode( wo.address_1 || ', Boston, MA ' || wo.zip_code, 1 ) AS g
    ON true;
  
-- VISUALIZE
SELECT *
FROM women_owned_geocoded;

-- SLIDE: K-nearest distance
-- Check which recent_building_permits is closest to each and draw a line between the two
SELECT wo.id AS wo_id, wo.business_name, 
    ST_MakeLine(wo.geog::geometry, rbo.geog::geometry) AS line, 
    rbo.*
FROM women_owned_geocoded AS wo 
    CROSS JOIN LATERAL 
    (SELECT rb.*,  rb.geog <-> wo.geog AS dist_m
        FROM recent_building_permits AS rb 
     ORDER BY dist_m LIMIT 1 ) AS rbo
ORDER BY rbo.geog <-> wo.geog; 

-- 3.0 ST_ASGeoJSON enhancement --
-- SLIDE: ST_AsGeoJSON can output a whole row including geometry
SELECT ST_AsGeoJSON(r)
FROM recent_building_permits AS r
LIMIT 3;


-- SLIDE: Build a feature collection from the features
SELECT json_build_object('type', 'FeatureCollection', 'features', 
    json_agg(ST_AsGeoJSON(r.*)::json) )
FROM recent_building_permits AS r;



--- WHAT IS NEW IN PostGIS 3.1 --

-- but first lets download boston neighborhoods:
-- https://data.boston.gov/dataset/boston-neighborhoods1/resource/fc2ef94d-b4fb-4780-8f4a-bcff7724e50a
DROP TABLE IF EXISTS boston_neighborhoods;
WITH d AS (SELECT  h.content::jsonb AS content
FROM http_get('https://bostonopendata-boston.opendata.arcgis.com/datasets/3525b0ee6e6b427f9aab5d0a1d0a1a28_0.geojson?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D') AS h)
SELECT
    r->'properties'->>'OBJECTID' AS object_id,
    r->'properties'->>'Name' AS name, 
    ST_GeomFromGeoJSON(r->>'geometry')::geography AS geog
INTO boston_neighborhoods
FROM d, jsonb_array_elements(d.content->'features') AS r;

SELECT *
FROM boston_neighborhoods
ORDER BY name;


-- SLIDE Grid building functions ST_HexagonGrid, ST_SquareGrid
-- 2249 is MA State Plane ft
-- hexagon of 0.005 degree edges
SELECT d.name, g.*
FROM boston_neighborhoods AS d, 
    ST_HexagonGrid(0.005,d.geog::geometry) AS g
WHERE ST_Intersects(d.geog, g.geom::geography);

-- square of 0.005 degree edges
SELECT d.name, g.*
FROM boston_neighborhoods AS d,
    ST_SquareGrid(0.005,d.geog::geometry) AS g
WHERE ST_Intersects(d.geog, g.geom::geography)
ORDER BY g.i, g.j, d.name;

-- SLIDE: Union / Intersection / ST_SubDivide and other functions now have precision
-- no precision 2249 is MA State Plane ft, 40 ms
SELECT ST_Union(ST_Transform(geog::geometry,2249))
FROM boston_neighborhoods;

-- Let's play with precision from 0 ft to 2000 ft in steps of 200 ft
SELECT i, ST_Union( ST_Transform(geog::geometry,2249), i ) AS geom_st_ft,
    ST_Transform(ST_Union( ST_Transform(geog::geometry,2249), i ),4326)::geography AS geog
FROM boston_neighborhoods, generate_series(0, 2000, 200) AS i
GROUP BY i;

-- SLIDE: Subdivide fixed precision
-- basic subdivide says max 200 vertices per geom below, 9 pieces
-- 2249 is MA State plane feet
SELECT sb.ord, sb.geom
FROM (SELECT ST_Transform(geog::geometry,2249) As geom,  
        ST_NPoints(geog::geometry) AS npoints
      FROM boston_neighborhoods WHERE name = 'Downtown') AS ref
      , ST_SubDivide(ref.geom,200) WITH ORDINALITY AS sb(geom,ord);
      
-- with 125ft precision requires GEOS 3.9
-- and max 200 vertices, note only 3 chunks
SELECT sb.ord, sb.geom
FROM (SELECT ST_Transform(geog::geometry,2249) As geom,  
        ST_NPoints(geog::geometry) AS npoints
      FROM boston_neighborhoods WHERE name = 'Downtown') AS ref
      , ST_SubDivide(ref.geom,200, 125) WITH ORDINALITY AS sb(geom,ord);

-- PostGIS 3.2 changes committed already
-- Several raster functions ST_InterpolateRaster, ST_Contour, ST_SetM, ST_SetZ, enhanced ST_Value

-- Build a raster using ST_InterpolateRaster
DROP TABLE IF EXISTS building_permits_raster;
CREATE TABLE building_permits_raster AS 
WITH p AS (
    -- theme points building permits
SELECT ST_Transform(geog::geometry, 2249) AS geom, 
        COUNT(*) AS num_permits , sum(total_fees) AS sum_fees 
     FROM recent_building_permits 
    GROUP BY geom
), 
ext AS (
    -- extent of our raster
    SELECT ST_SetSRID(ST_Extent(geom),2249) AS geom
    FROM p
)
SELECT -- https://postgis.net/docs/RT_ST_InterpolateRaster.html
      ST_InterpolateRaster(
        -- collect the points and fees as the z-value
        ST_Collect( 
          ST_Force3D(
              p.geom ,
              sum_fees
              )
        ),
   'linear',
    ST_AddBand( 
        ST_MakeEmptyRaster(width => ((ST_XMax(ext.geom) - ST_XMin(ext.geom))/50)::integer,
                           height => ((ST_YMax(ext.geom) - ST_YMin(ext.geom))/50)::integer, 
                           upperleftx => ST_XMin(ext.geom),
                           upperlefty => ST_YMax(ext.geom),
                           scalex => 50,
                           scaley => - 50,
                           skewx => 0, skewy =>0, 
                           srid => 2249)
        , '16BSI' )
  ) AS rast
FROM p, ext
GROUP BY ext.geom;

SELECT ST_Value(rast, ST_Transform(ST_SetSRID(ST_Point(-71.0573, 42.3583), 4326),2249)  )
FROM building_permits_raster;

-- spot check to make sure it covers the area you think it does
SELECT ST_Transform(rast::geometry,4326)
FROM building_permits_raster;

-- check how much cost of approximate building permits in my area
-- gives 1315
SELECT ST_Value(rast
                  , ST_Transform(ST_SetSRID(ST_Point(-71.0573, 42.3583), 4326),2249)
               )
FROM building_permits_raster;

-- spot check from original permit data, not different because raster is linear disperse
SELECT AVG(total_fees)
FROM recent_building_permits where ST_DWithin(geog,ST_Point(-71.0573, 42.3583)::geography,100);


-- output as a png, base64 (copy and paste contents into text file and view in web browser)
SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
WITH boston AS (SELECT ST_Union( ST_Transform(geog::geometry,2249)  ) AS geom
               FROM boston_neighborhoods
               )
SELECT '<img src="data:image/png;base64,' 
    || encode(ST_AsPNG( 
    -- https://postgis.net/docs/RT_ST_ColorMap.html
        ST_ColorMap(
        ST_Clip(rast, boston.geom ) 
        , '100% 255 0 0 255
90% 200 0 5 255
50% 100 0 50 255
25% 50 0 220 255
10% 2 0 250 255
1% 0 0 255 255
nv 255 255 255 255'
    
    ) ), 'base64') || '" />'  
FROM building_permits_raster, boston;