/**
Extensions to cover:
postgis
http (not a spatial extension, but works great with everything)
postgis_sfcgal
postgis_raster
ogr_fdw
address_standardizer
postgis_tiger_geocoder
h3, h3_postgis
pointcloud, pointcloud_postgis

**/
-- https://postgis.us/presentations/postgisday_2023.sql
SELECT version();


--https://postgis.us/presentations/postgisday2023.html
CREATE DATABASE gisdb;
ALTER DATABASE gisdb SET search_path=postgisday2023,contrib,postgis,tiger,topology,public;
\c gisdb
CREATE SCHEMA IF NOT EXISTS contrib;
DROP SCHEMA IF EXISTS postgisday2023 CASCADE;
CREATE SCHEMA IF NOT EXISTS postgisday2023;

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

-- In PostGIS 3.3, some guy brought us ST_Letters, https://postgis.net/docs/ST_Letters.html
-- allowing us to say hello in PostGIS and even see hello on a map
SELECT ST_SetSRID(ST_Letters('hello'),4326);

CREATE EXTENSION IF NOT EXISTS postgis_sfcgal;

SELECT sk.geom, h.geom, ST_Collect(sk.geom, h.geom)
FROM ST_Letters('H') AS h(geom), ST_StraightSkeleton(h.geom ) AS sk(geom);

-- but immediately the crowd wanted more. They wanted to use their own fonts.
-- They wanted to say hello in different languages.
-- BUT I'll accept another argument to this function, and he ended with these famous fighting words
-- "Figure out what to do with it yourself"

-- and so I bring you one way you can create your own font set with the power of postgis and postgis_raster 

-- START: Make a star wars font set using http, postgis_raster, and postgis
-- First we start off with an extension created by that same guy
-- you may need to compile your own - https://github.com/pramsey/pgsql-http
-- not a spatial extension, but tastes great with lots of other extensions
CREATE EXTENSION IF NOT EXISTS http SCHEMA contrib;
CREATE SCHEMA IF NOT EXISTS postgis;
-- will install both postgis and postgis_raster in schema postgis
CREATE EXTENSION IF NOT EXISTS postgis_raster SCHEMA postgis CASCADE;

-- https://postgis.net/docs/en/postgis_gdal_enabled_drivers.html
SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';

-- for future connections
ALTER DATABASE gisdb SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
-- check drivers for read / write
SELECT * FROM ST_GDALDrivers() ORDER BY long_name;

DROP TABLE IF EXISTS raster_raw;
CREATE TEMP TABLE raster_raw(name text, rast raster);
-- if get error ssl cert
--SELECT http_set_curlopt('CURLOPT_CAINFO', 'C:/Program Files/PostgreSQL/16/ssl/certs/ca-bundle.crt');
INSERT INTO raster_raw(name, rast)
SELECT 'star wars', ST_FromGDALRaster(text_to_bytea(h.content))
FROM http_get('https://i.etsystatic.com/25651615/r/il/0a1e81/3902154889/il_794xN.3902154889_cx7x.jpg') AS h;

-- https://as2.ftcdn.net/v2/jpg/04/31/23/73/1000_F_431237393_jyPB10KLUlKdtgmI3xLF5eh7zQhB5g9f.jpg

-- view in browser
SELECT 'data:image/png;base64,' || encode(ST_AsPNG(
    ST_Band(rast,1) ),'base64') 
  FROM raster_raw
	WHERE name = 'star wars' ;

-- check range of pixel values, band 1 into 3 buckets
SELECT hg.*
FROM raster_raw , ST_Histogram(rast,1, 3) AS hg
WHERE name = 'star wars';

-- reclassify into black and white https://postgis.net/docs/en/RT_ST_Reclass.html
INSERT INTO raster_raw(name, rast)
SELECT 'star wars reclass', ST_Reclass(ST_Band(rast,1), 1, '44-184: 1, 185-:0', '1BB', 0)
FROM raster_raw
WHERE name = 'star wars';

-- check width and height, original 794, 650
SELECT name, ST_Width(rast), ST_Height(rast)
FROM raster_raw;

-- clip to minimum nodata boundary https://postgis.net/docs/en/RT_ST_MinConvexHull.html
UPDATE raster_raw SET rast = ST_Clip(ST_Band(rast,1), ST_MinConvexHull(rast));

-- check width and height, new 700, 516
SELECT name, ST_Width(rast), ST_Height(rast)
FROM raster_raw;

-- check reclass as polygon
SELECT ST_Polygon(rast)
FROM raster_raw WHERE name = 'star wars reclass';


-- we want each tile to be 500 x 500 and there are 4 across and 3 down
UPDATE raster_raw SET rast = ST_Resize(rast, 500*4,500*3, algorithm => 'cubic'  )
WHERE name = 'star wars reclass';
	
DROP TABLE IF EXISTS stencil_set; 
CREATE TABLE stencil_set(name text, geom geometry, rast raster);
-- m starts at 109, and m is hardcoded as the spacer test so we need an m in our fontset
-- https://postgis.net/docs/en/RT_ST_Tile.html
INSERT INTO stencil_set(name, geom, rast)
-- use CHR (built in postgresql function) to convert numbers to letters						
SELECT CHR(t.ordinality::integer + 100), f.geom, t.rast
-- use ST_Tile to chop our stencil into 500x500 pieces
FROM raster_raw, ST_Tile(rast,1,500,500 ) WITH ORDINALITY AS t(rast), ST_Polygon(t.rast ) AS f(geom)
WHERE name = 'star wars reclass' ;

SELECT name, geom
FROM stencil_set;


-- get rid of boxes and other junk and straigten out the jagged edges
-- one approach
SELECT name, ST_SimplifyVW(ST_Union(ST_BuildArea(ST_ExteriorRing(dp.geom))),5), MAX(ST_Area(dp.geom))
FROM stencil_set, ST_Dump(geom) AS dp
WHERE ST_Area(dp.geom) BETWEEN 3 AND 4000 or ST_Area(dp.geom) > 6000 
GROUP BY name
ORDER BY name;

-- this approach is more visually appealing
-- get rid of boxes and other junk and straigten out the jagged edges
-- https://postgis.net/docs/en/ST_SimplifyPolygonHull.html (sorry only GEOS 3.11+ has this)
SELECT name, ST_Union(ST_BuildArea(ST_ExteriorRing(dp.geom))) As original, ST_SimplifyPolygonHull( ST_Union(ST_BuildArea(ST_ExteriorRing(dp.geom)))
	,vertex_fraction => 0.28, is_outer => false), MAX(ST_Area(dp.geom))
FROM stencil_set, ST_Dump(geom) AS dp
WHERE ST_Area(dp.geom)  BETWEEN 3 AND 4050 or ST_Area(dp.geom) > 6000 
GROUP BY name
ORDER BY name;

-- okay lets back update, once we've fiddled with area and simplification
WITH a AS (
SELECT name, ST_Union(ST_BuildArea(ST_ExteriorRing(dp.geom))) As original, ST_SimplifyPolygonHull( ST_Union(ST_BuildArea(ST_ExteriorRing(dp.geom)))
	,vertex_fraction => 0.28, is_outer => false) AS geom, MAX(ST_Area(dp.geom))
FROM stencil_set, ST_Dump(geom) AS dp
WHERE ST_Area(dp.geom)  BETWEEN 3 AND 4050 or ST_Area(dp.geom) > 6000 
GROUP BY name
ORDER BY name
)
UPDATE stencil_set AS n SET geom = a.geom
	FROM a
WHERE a.name = n.name;


-- build a font-set from the stencil
DROP TABLE IF EXISTS font_set;
CREATE TABLE font_set(name text, font json);
INSERT INTO font_set(name,font)
SELECT 'star wars', json_object(array_agg(name),
-- https://postgis.net/docs/en/ST_AsTWKB.html The tiny well-known binary, a light weight fantastic format you probably didn't know existed,
-- and that ST_Letters needs for fonts
-- translate everything to start at 0,0, otherwise goofy stuff happens when y is different, cause ST_Letters has bugs
-- also you need an m, cause the guy couldn't imagine a world without the letter m, that's why we had all that nonsense with CHR
	array_agg(replace(encode(ST_AsTWKB( ST_Translate(geom, -ST_XMax(geom), -ST_YMax(geom))),'base64'), E'\n', '')) )
FROM stencil_set;

-- Imagine what masterful secret messages we could pass to each other, 
-- if only we understood more than kindergarden cryptography
SELECT f.geom
FROM ST_Letters('hello', (SELECT font FROM font_set WHERE name = 'star wars')) AS f(geom);

-- END: Make a star wars font set using http, postgis_raster, and postgis


-- START: ogr_fdw https://github.com/pramsey/pgsql-ogr-fdw
CREATE EXTENSION IF NOT EXISTS ogr_fdw SCHEMA postgis;
DROP SCHEMA  IF EXISTS staging CASCADE;
CREATE SCHEMA staging;


-- check drivers you have:
SELECT name FROM unnest(ogr_fdw_drivers()) AS f(name) ORDER BY name;
 
 
			
-- original location is https://data.humdata.org/dataset/9159a100-8672-4c5b-8162-c24459514bfa/resource/3f4ee676-1936-45e5-9a6a-aa224cf3aa39/download/kontur_population_20220630_r4.gpkg.gz (but redirects make it not work)
-- details about the dataset here - https://www.kontur.io/portfolio/population-dataset/
-- NOTE: We are reading a gzipped file of a geopackage db, across the internet, how cool is that?
 CREATE SERVER fds_humapop
    FOREIGN DATA WRAPPER ogr_fdw
    OPTIONS(
        datasource  '/vsigzip//vsicurl/https://postgis.us/downloads/postgisday2023/kontur_population_20220630_r4.gpkg.gz' 
            , format 'GPKG');
	
CREATE SCHEMA IF NOT EXISTS staging;	
-- bring in all the tables (there is only one)
IMPORT FOREIGN SCHEMA ogr_all FROM SERVER fds_humapop INTO staging;	
-- note we now have a population table -- hold that thought
SELECT table_name
 FROM information_schema.tables 
 WHERE table_schema = 'staging';
 
-- END: ogr_fdw

-- START: Address standardizer (new in 3.4, debug_address_standizer)
CREATE EXTENSION IF NOT EXISTS address_standardizer SCHEMA contrib;
CREATE EXTENSION IF NOT EXISTS address_standardizer_data_us SCHEMA contrib;

-- https://postgis.net/docs/standardize_address.html
SELECT *
FROM standardize_address('us_lex',
                'us_gaz', 'us_rules', '51 Franklin Street, Fifth Floor', 'Boston, MA 02110-1301, USA');
				
-- START: debug the rules https://postgis.net/docs/debug_standardize_address.html
-- look at the raw json new in PostGIS 3.4.0
SELECT ds::json
FROM debug_standardize_address('us_lex',
                'us_gaz', 'us_rules', '51 Franklin Street, Fifth Floor', 'Boston, MA 02110-1301, USA') AS ds;
	
-- look at the rules
SELECT itb As raw_jsonb, itb->>'no', itb->>'score' AS score
FROM jsonb(debug_standardize_address('us_lex',
                'us_gaz', 'us_rules', '51 Franklin Street, Fifth Floor', 'Boston, MA 02110-1301, USA')) AS s
				CROSS JOIN jsonb_array_elements(s->'rules') AS it, jsonb(it) AS itb;
				
-- each rule
-- look at the input tokens
SELECT ds::json
FROM debug_standardize_address('us_lex',
                'us_gaz', 'us_rules', '51 Franklin Street, Fifth Floor', 'Boston, MA 02110-1301, USA') AS ds;
SELECT it->>'pos' AS position, it->>'word' AS word, it->>'stdword' AS standardized_word,
            it->>'token' AS token, it->>'token-code' AS token_code
FROM jsonb(debug_standardize_address('us_lex',
                'us_gaz', 'us_rules', '51 Franklin Street, Fifth Floor', 'Boston, MA 02110-1301, USA')) AS s
				CROSS JOIN jsonb_array_elements(s->'input_tokens') AS it;

-- END: Address standardizer (new in 3.4, debug_address_standizer)

-- START tiger_geocoder (switch to database that has data loaded)
CREATE EXTENSION IF NOT EXISTS postgis_tiger_geocoder;
-- load data following these instructions - https://postgis.net/docs/manual-3.4/en/postgis_installation.html#install_tiger_geocoder_extension
-- profit
SELECT ST_SetSRID(geomout,4326), ST_AsText(geomout)
from geocode('51 Franklin Street, Fifth Floor, Boston, MA 02110-1301', 1);


SELECT ST_SetSRID(r.intpt[1],4326) , addy[1].*, street[1], street[2], *
FROM reverse_geocode(ST_Point(-71.05871, 42.35555)) AS r;

SELECT g.*, ST_SetSRID(geomout,4326)
FROM geocode_intersection('Hawley St', 'Franklin St','MA', '02110') AS g;

-- END: tiger_geocoder

-- START: make uber hexagons with h3, h3_postgis h3-pg https://github.com/zachasme/h3-pg
-- the CASCADE forces install of h3 extension as well
CREATE EXTENSION IF NOT EXISTS h3_postgis SCHEMA postgis CASCADE;

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

-- remember that population data, it's slow cause it's not in our database
-- note the h3 column is just text, but it's canonical textural representation of an h3index
SELECT * FROM staging.population limit 1;

-- now it's in our database
-- 70296 rows
CREATE TABLE human_pop AS 
SELECT h3::h3index, population
FROM staging.population;

CREATE INDEX ix_human_pop_h3 ON human_pop USING brin(h3);

-- what is the population near me?
-- what's the population around here? Lets have h3 answer this question
-- this is a very low res pop set, pick a different h3 and you can zoom in more and aggregate
SELECT  h3_cell_to_boundary_geometry(p.h3), p.population
FROM human_pop AS p
WHERE p.h3 && h3_lat_lng_to_cell(ST_Point(-71.05871,42.35555,4326), 10);

-- lets try the new raster support
SELECT
    h3_cell_to_boundary_geometry(hrs.h3),
	(h3_raster_summary_stats_agg(hrs.stats)).*
FROM
    stencil_set AS r, ST_Resize(r.rast,0.5,0.5) AS rs(rast),
    h3_raster_summary(ST_SetSRID(ST_SetUpperLeft(rs.rast,0,ST_Height(rs.rast)) ,4326),3) AS hrs
WHERE r.name = 'e' 
GROUP BY hrs.h3;

-- END: Make uber hexagons with h3, h3_postgis

-- 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) -- https://pdal.io/
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
-- stop here