-- https://postgis.us/presentations/postgisday_2022.sql
SELECT version();
--https://postgis.us/presentations/postgisday2022.html
CREATE DATABASE gisdb;
ALTER DATABASE gisdb SET search_path=postgisday2022,contrib,postgis,tiger,topology,public;
\c gisdb
CREATE SCHEMA IF NOT EXISTS contrib ;
--DROP SCHEMA IF EXISTS postgisday2022 CASCADE;
CREATE SCHEMA IF NOT EXISTS postgisday2022;

SET search_path=postgisday2022,contrib,postgis,tiger,topology,public;

-- Install postgis 
CREATE SCHEMA IF NOT EXISTS postgis;
CREATE EXTENSION IF NOT EXISTS postgis SCHEMA postgis;

SELECT postgis_extensions_upgrade();
SELECT postgis_full_version();

-- New in PostGIS 3.3.
-- https://postgis.net/docs//manual-3.3/ST_Letters.html
-- geometry ST_Letters(text letters, optional json font);
-- hello world in flat earth
SELECT ST_SetSRID(ST_letters('Hello World'),4326);

-- hello world in spheroidal earth, notice the wrap around
SELECT ST_SetSRID(ST_letters('Hello World'),4326)::geography;


-- accepts fontsets with letters encoded in TWKB (tiny twkb format) and key value in json
-- fontset borrowed from Bruce Rindahl and reformatted to twkb
SELECT fs.font_name,  ST_SetSRID(ST_Letters('Hello World', fs.letters),4326)
    FROM letters_font_set AS fs;


-- address standardizer
CREATE EXTENSION address_standardizer SCHEMA postgis;

-- US data set for standardizing
CREATE EXTENSION IF NOT EXISTS address_standardizer_data_us SCHEMA postgis;
-- uses gazetteer, lexicon, and rules - https://postgis.net/docs/manual-3.2/standardize_address.html

SELECT * 
FROM json_each_text( to_json ( standardize_address(
    'us_lex',
    'us_gaz', 
    'us_rules', 
    '51 Franklin Street, Fifth Floor, Boston, MA 02110') ) );

-- install geocoder -- prep data https://postgis.net/docs/postgis_installation.html#install_tiger_geocoder_extension
CREATE EXTENSION postgis_tiger_geocoder CASCADE;


SELECT pprint_addy(g.addy), ST_SetSRID(g.geomout, 4326), ST_AsText(g.geomout), g.rating 
 FROM geocode('51 Franklin Street, Fifth Floor, Boston, MA 02110',1) AS g;


DROP TABLE IF EXISTS test_geometries;
CREATE TABLE test_geometries(id integer GENERATED BY DEFAULT AS IDENTITY, guid uuid NOT NULL DEFAULT gen_random_uuid ()  PRIMARY KEY,
                                name text, geom geometry(geometry, 4326) );
                                
INSERT INTO test_geometries(name, geom)
SELECT 'Hello Boston',  l.geom
FROM geocode('51 Franklin Street, Fifth Floor, Boston, MA 02110', 1) AS gc, 
    ST_Translate(ST_SetSRID(ST_Scale(ST_Letters('Hello World'),0.001,0.001), 4326), ST_X(gc.geomout), ST_Y(gc.geomout) )  AS l(geom);
    

INSERT INTO test_geometries(name, geom)
SELECT 'Boston Again', ST_SetSRID(the_geom,4326)
FROM tiger.place
WHERE name = 'Boston';
    
SELECT *, ST_AsText(geom) FROM test_geometries;

-- h3-pg https://github.com/zachasme/h3-pg
-- binaries for windows for PG 15 -- https://winnie.postgis.net/download/windows/pg15/buildbot/extras/h3pg-pg15-binaries-mainw64.zip
CREATE EXTENSION IF NOT EXISTS h3 SCHEMA postgis;
CREATE EXTENSION IF NOT EXISTS h3_postgis SCHEMA postgis;


-- check version you are running
SELECT h3_get_extension_version();

-- add a point index column
ALTER TABLE test_geometries ADD h3 h3index GENERATED ALWAYS AS ( h3_lat_lng_to_cell( ST_Centroid(geom), 5) ) STORED;

-- convert a h3 to a postgis point
SELECT ST_Collect(geom, h3::geometry) FROM test_geometries;

DROP TABLE IF EXISTS test_h3_index;
-- 632 ms 3735 rows
CREATE TABLE test_h3_index(guid uuid, h3_key h3index, CONSTRAINT pk_h3_index PRIMARY KEY (h3_key, guid) );
INSERT INTO test_h3_index(guid, h3_key)
SELECT DISTINCT ON(g.id, h3) g.guid, h3.h3
FROM test_geometries AS g 
    CROSS JOIN h3_polygon_to_cells(g.geom,9) AS h3;
    
SELECT guid, h3_cell_area(h3_key, 'rads^2') AS area_rads,  h3_cell_area(h3_key, 'm^2') AS area_meters
FROM test_h3_index;

SELECT h3_get_hexagon_edge_length_avg(9,'m'); -- 174.3 meters
SELECT h3_cell_area('892a3076bcfffff'::h3index, 'm^2');
    
CREATE INDEX ix_test_h3_index_h3_key ON test_h3_index USING brin(h3_key);

-- h3_cell_to_boundary_geometry returns hexagon, casting to geometry returns the centroid of hexagon
SELECT ST_Collect(h3_cell_to_boundary_geometry(h3_key)), COUNT(1), ST_Collect(h3_key::geometry)
FROM test_h3_index AS h3 INNER JOIN test_geometries AS g ON h3.guid = g.guid
WHERE g.id = 2;


SELECT 174.3/110574; --0.0027631993054425091
DROP TABLE IF EXISTS test_postgis_hexagons;
-- 3725 rows, 181 ms
CREATE TABLE test_postgis_hexagons(guid uuid,  i integer, j integer, geom geometry(POLYGON,4326), PRIMARY KEY (i, j, guid));
INSERT INTO test_postgis_hexagons(guid, i, j, geom)
SELECT g.guid, h.i, h.j, h.geom
FROM test_geometries AS g INNER JOIN ST_HexagonGrid(0.00235,g.geom) AS h ON ST_Intersects(g.geom, h.geom);

SELECT ST_Collect(h.geom), COUNT(1), ST_Collect(ST_Centroid(h.geom))
FROM test_postgis_hexagons AS h INNER JOIN test_geometries AS g ON h.guid = g.guid
WHERE g.id = 2;

-- native postgresql polygon type
WITH a AS (SELECT (h3_cells_to_multi_polygon(array_agg(h3_key))).*
FROM test_h3_index
GROUP BY guid)
SELECT exterior, holes, exterior::geometry --postgis has cast for native polygon types
FROM a;


-- intersect
SELECT SUM(h3_cell_area(i1.h3_key, 'm^2')) AS intersecting_area,  ST_Collect(h3_cell_to_boundary_geometry(i1.h3_key)) AS hexagons
FROM test_h3_index AS i1 INNER JOIN test_h3_index AS i2 ON ( i1.h3_key = i2.h3_key )
WHERE i1.guid <> i2.guid


-- convert back to polygon
SELECT h3_cells_to_multi_polygon_geometry(array_agg(h3_key)) --, ST_Collect(h3_key::geometry)
FROM test_h3_index
GROUP BY guid;
    
-- ogr_fdw source at: https://github.com/pramsey/pgsql-ogr-fdw 
CREATE EXTENSION IF NOT EXISTS ogr_fdw SCHEMA postgis;
CREATE SCHEMA IF NOT EXISTS staging ;

-- see what drivers you have
SELECT driver FROM  unnest(ogr_fdw_drivers()) AS driver ORDER BY driver;

DROP SERVER IF EXISTS svr_shps CASCADE;
-- set to false cause there's a bug that it tries to see if it is readable https://github.com/pramsey/pgsql-ogr-fdw/issues/221
CREATE SERVER IF NOT EXISTS svr_shps 
 FOREIGN DATA WRAPPER ogr_fdw
  OPTIONS (datasource '/vsizip//vsicurl/http://postgis.us/downloads/postgisday2022/data.zip',format 'ESRI Shapefile',
          updateable 'false'
          );

IMPORT FOREIGN SCHEMA ogr_all FROM SERVER svr_shps INTO staging;
SELECT *
FROM staging.boston_neighborhoods;

DROP SERVER IF EXISTS fds_gpx CASCADE;
CREATE SCHEMA IF NOT EXISTS staging; 

-- trace from here https://www.openstreetmap.org/trace/5360016/data
CREATE SERVER fds_gpx
    FOREIGN DATA WRAPPER ogr_fdw
    OPTIONS(
        datasource  '/vsicurl/http://postgis.us/downloads/postgisday2022/5360016.gpx' 
            , format 'GPX');
-- ogr_all is catchall for give me everything in the data source
IMPORT FOREIGN SCHEMA ogr_all FROM SERVER fds_gpx INTO staging;
SELECT table_name
 FROM information_schema.tables 
 WHERE table_schema = 'staging';
 
-- mobilitydb 
-- https://mobilitydb.com/
SHOW shared_preload_libraries;
-- ALTER SYSTEM SET shared_preload_libraries = 'postgis-3';
-- need to restart postgresql service for above to take effect
DROP EXTENSION mobilitydb CASCADE; --this destroys data so be careful
CREATE EXTENSION IF NOT EXISTS mobilitydb SCHEMA postgis;
-- show what you are running -- MobilityDB 1.1.0, PostgreSQL 15.1 (Ubuntu 15.1-1.pgdg22.04+1), PostGIS 3.3.0
SELECT mobilitydb_full_version();

DROP TABLE IF EXISTS gpx_trips;
CREATE TABLE gpx_trips(file_name text, 
                       trip tgeogpoint, 
                       CONSTRAINT  pk_gpx_trips PRIMARY KEY(file_name)  );
                       
-- convert each geom time pair into an instant then aggregate into a geog point                       
INSERT INTO gpx_trips(file_name, trip)
SELECT '5360016.gpx',   tgeogpoint_seq( array_agg( sdp.instant ORDER BY sdp.ts))
FROM 
    (SELECT DISTINCT ON(time) tgeogpoint_inst(geom::geography, (time::text || ' UTC')::timestamptz ) AS instant, time AS ts
     FROM 
         staging.track_points
     ORDER BY time, geom
    ) AS sdp;

CREATE INDEX ix_gpx_trips_trip_gist ON gpx_trips USING gist(trip);
-- cast to a PostGIS LINESTRING M (linestring with measure in epoch) SELECT to_timestamp(1653343114)
SELECT file_name, trip, trip::geography AS path, duration(trip), 
    speed(trip) AS speed, 
    ST_AsText(trip::geography) AS mgpath
FROM gpx_trips
WHERE trip && '[2022-11-15, 2022-11-16)'::period;

-- aggregate functions
SELECT file_name, 
    twavg(speed(trip)) AS weigted_average_speed,
    tcount(trip)
FROM gpx_trips 
GROUP BY file_name;

-- create a table for our worlds
DROP TABLE IF EXISTS postgisday2022.hw;
CREATE TABLE postgisday2022.hw(name text PRIMARY KEY, geom geometry(MULTIPOLYGON,4326));

INSERT INTO postgisday2022.hw(name, geom)
WITH a AS (
    SELECT ST_Scale(
       -- hello world
        ST_Letters('Hello World'), 
       -- scaling to 10%
        ST_Point(0.1,0.1) 
    ) AS geom
   )
SELECT 'Hello World' AS name,  ST_SetSRID(
       a.geom
    -- on earth wgs 84 long lat
    , 4326)
    -- enforce multipolygon earth
    ::geometry(MULTIPOLYGON, 4326) AS geom
FROM a;

-- convert a whole row to json
SELECT ST_AsGeoJSON(hw)::jsonb
FROM postgisday2022.hw AS hw;


-- # postgis topology of Hello world
DROP EXTENSION IF EXISTS postgis_topology CASCADE;
CREATE EXTENSION IF NOT EXISTS postgis_topology CASCADE;
DROP SCHEMA IF EXISTS postgisday2022_topo CASCADE;

-- create a new topology
-- https://postgis.net/docs/CreateTopology.html
-- integer CreateTopology(varchar topology_schema_name, integer srid, double precision prec);
SELECT topology.CreateTopology('postgisday2022_topo', 4326, 0.0001);

-- https://postgis.net/docs/AddTopoGeometryColumn.html
-- integer AddTopoGeometryColumn(varchar topology_name, varchar schema_name, varchar table_name, varchar column_name, varchar feature_type);
-- returns layer id of new topology
SELECT topology.AddTopoGeometryColumn(
	'postgisday2022_topo',
	'postgisday2022', 'hw', 'topo', 'POLYGON'); 
    
UPDATE postgisday2022.hw SET topo = topology.toTopoGeom(geom, 'postgisday2022_topo',1);

SELECT topo, (topo).*, topo::geometry
FROM postgisday2022.hw;

-- lets check our edges
SELECT * FROM postgisday2022_topo.edge;

-- START PGROUTING
-- # pgRouting in words
-- https://docs.pgrouting.org/
DROP EXTENSION IF EXISTS pgrouting;
CREATE EXTENSION IF NOT EXISTS pgrouting SCHEMA postgis;
-- check version you are running
SELECT pgr_version();
-- full details of what dependencies you are running with
SELECT *
FROM pgr_full_version();

-- https://docs.pgrouting.org/latest/en/pgr_createTopology.html

-- lets make a path
DROP TABLE IF EXISTS pgr_roads2;
DROP TABLE IF EXISTS pgr_roads2_vertices_pgr;
CREATE TABLE pgr_roads2(id bigint GENERATED ALWAYS AS IDENTITY, name text, geom geometry(LINESTRING,4326), source bigint, target bigint, 
                        cost float,x1 float, y1 float, x2 float, y2 float);
-- create dummy data - 70 roads
INSERT INTO pgr_roads2(name, geom, cost)
SELECT 'road ' || dp.path[1] AS name, ST_SnaptoGrid(dp.geom,0.00001) AS geom, ST_Length(dp.geom::geography)
FROM 
    ST_SetSRID( 
            ST_Translate(
                    ST_Scale( 
                            ST_Letters('Hello World'),0.25, 0.25), -71.05871, 42.3555)
               , 4326) AS hw(geom)  
CROSS JOIN ST_Boundary(hw.geom) AS g(g)
CROSS JOIN ST_Union(g.g, ST_LongestLine(g.g,g.g)) AS gp
CROSS JOIN ST_Dump(gp) AS dp;


-- https://docs.pgrouting.org/latest/en/pgr_createTopology.html
-- https://docs.pgrouting.org/latest/en/pgr_createVerticesTable.html
-- https://docs.pgrouting.org/latest/en/pgr_extractVertices.html#id12 #newer and faster
CREATE TABLE pgr_roads2_vertices_pgr AS
SELECT  * FROM pgr_extractVertices(
  'SELECT id, geom FROM pgr_roads2');
-- update source  
WITH
out_going AS (
  SELECT id AS vid, unnest(out_edges) AS eid, x, y
  FROM pgr_roads2_vertices_pgr
)
UPDATE pgr_roads2
SET source = vid, x1 = x, y1 = y
FROM out_going WHERE id = eid;

-- update target 
WITH
in_coming AS (
  SELECT id AS vid, unnest(in_edges) AS eid, x, y
  FROM pgr_roads2_vertices_pgr
)
UPDATE pgr_roads2
SET target = vid, x2 = x, y2 = y
FROM in_coming WHERE id = eid;

-- https://github.com/pgrouting (available in yum.postgresql.org, apt.postgresql.org, and windows postgis bundle)
-- https://docs.pgrouting.org/3.3/en/pgr_connectedComponents.html
-- pgr_connectedComponents(edges_sql)
-- lists all components reachable from each other
WITH cc AS (SELECT cc.*, n.geom 
FROM pgr_connectedComponents(
    'SELECT id, source, target, cost FROM pgr_roads2'
) AS cc INNER JOIN pgr_roads2_vertices_pgr AS n ON cc.node = n.id
)
SELECT cc.component, 
    ST_Collect( cc.geom ORDER BY cc.seq) AS geom,
    array_agg(DISTINCT cc.seq ORDER BY cc.seq) AS nodes
FROM cc 
GROUP BY cc.component
UNION ALL
-- overlay the roads
SELECT -1, ST_Collect(geom), null
FROM pgr_roads2
ORDER BY 1;
-- https://docs.pgrouting.org/3.3/en/pgr_dijkstra.html
-- pgr_dijkstra(Edges SQL, start_vid,  end_vid  [, directed])
-- route 
WITH path AS (SELECT d.* , n.geom AS node_geom, e.geom As edge_geom
FROM pgr_dijkstra(
    'SELECT id, source, target, cost FROM pgr_roads2',
   1, 18, directed => false
) AS d LEFT JOIN pgr_roads2 AS e ON d.edge = e.id
    LEFT JOIN pgr_roads2_vertices_pgr AS n ON d.node = n.id 
   ORDER BY d.seq    )
SELECT 0, ST_Union(edge_geom)
FROM path
UNION ALL
SELECT seq, node_geom
FROM path;

-- END PGROUTING

-- START PGPOINTCLOUD
-- pgpointcloud https://github.com/pgpointcloud/  (available in yum.postgresql.org, apt.postgresql.org, and windows postgis bundle)
-- docs: https://pgpointcloud.github.io/pointcloud/
CREATE EXTENSION IF NOT EXISTS pointcloud SCHEMA postgis;
-- contains functions for casting between pcpatches and postgis geometries
CREATE EXTENSION IF NOT EXISTS pointcloud_postgis SCHEMA postgis;

-- use pdal to load real data - https://www.pdal.io/

--see formats
SELECT * FROM pointcloud_formats;

-- add a format (this is normally done by pdal)
INSERT INTO pointcloud_formats (pcid, srid, schema) VALUES (1, 4326,
'<?xml version="1.0" encoding="UTF-8"?>
<pc:PointCloudSchema xmlns:pc="http://pointcloud.org/schemas/PC/1.1"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <pc:dimension>
    <pc:position>1</pc:position>
    <pc:size>4</pc:size>
    <pc:description>X coordinate as a long integer. You must use the
                    scale and offset information of the header to
                    determine the double value.</pc:description>
    <pc:name>X</pc:name>
    <pc:interpretation>int32_t</pc:interpretation>
    <pc:scale>0.01</pc:scale>
  </pc:dimension>
  <pc:dimension>
    <pc:position>2</pc:position>
    <pc:size>4</pc:size>
    <pc:description>Y coordinate as a long integer. You must use the
                    scale and offset information of the header to
                    determine the double value.</pc:description>
    <pc:name>Y</pc:name>
    <pc:interpretation>int32_t</pc:interpretation>
    <pc:scale>0.01</pc:scale>
  </pc:dimension>
  <pc:dimension>
    <pc:position>3</pc:position>
    <pc:size>4</pc:size>
    <pc:description>Z coordinate as a long integer. You must use the
                    scale and offset information of the header to
                    determine the double value.</pc:description>
    <pc:name>Z</pc:name>
    <pc:interpretation>int32_t</pc:interpretation>
    <pc:scale>0.01</pc:scale>
  </pc:dimension>
  <pc:dimension>
    <pc:position>4</pc:position>
    <pc:size>2</pc:size>
    <pc:description>The intensity value is the integer representation
                    of the pulse return magnitude. This value is optional
                    and system specific. However, it should always be
                    included if available.</pc:description>
    <pc:name>Intensity</pc:name>
    <pc:interpretation>uint16_t</pc:interpretation>
    <pc:scale>1</pc:scale>
  </pc:dimension>
  <pc:metadata>
    <Metadata name="compression">dimensional</Metadata>
  </pc:metadata>
</pc:PointCloudSchema>');

DROP TABLE IF EXISTS hw_pcpoints;
CREATE TABLE hw_pcpoints (
    id bigint GENERATED ALWAYS AS IDENTITY PRIMARY KEY,
    pt PCPOINT(1) -- the one means use format with id = 1
);

-- A table of patches
DROP TABLE IF EXISTS hw_pcpatches;
CREATE TABLE hw_pcpatches (
    id bigint GENERATED ALWAYS AS IDENTITY PRIMARY KEY,
    pa PCPATCH(1) -- the one means use format with id = 1
);

SELECT * FROM pointcloud_columns;
-- https://postgis.net/docs/ST_GeneratePoints.html
INSERT INTO hw_pcpoints(pt)
SELECT PC_MakePoint(1, ARRAY[ST_X(gd.geom), ST_Y(gd.geom), random()*1000, random()*40])
FROM ST_Dump(
        ST_GeneratePoints( 
            ST_Scale(ST_Letters('Hello World'), 0.2, 0.2), 1000, 1)
    ) AS gd;
    
-- a patch is an array of pcpoints
-- https://pgpointcloud.github.io/pointcloud/functions/patchs.html#pc-patch
INSERT INTO hw_pcpatches(pa)
SELECT PC_Patch(pt)
FROM hw_pcpoints
GROUP BY id/100;

SELECT ST_Force2D(PC_EnvelopeGeometry(pa))
FROM hw_pcpatches;

CREATE INDEX ON hw_pcpatches USING gist(PC_EnvelopeGeometry(pa) gist_geometry_ops_nd);
SELECT hp.id, PC_AsText(subpatch), PC_PatchAvg(hp.pa, 'Intensity'), subpatch::geometry
FROM hw_pcpatches AS hp , PC_FilterEquals(hp.pa, 'Intensity', 1) AS subpatch;
-- END PGPOINTCLOUD

