Skip to content

Home

Jump to bottom
Frank Horowitz edited this page May 3, 2019 · 13 revisions

Note added later: I finally found a decent book that describes building SQL queries in comprehensible language. Have a look at SQL Queries For Mere Mortals. (Unfortunately, I was impatient and bought the Kindle edition, so I have no physical book to lend others.) But after spending some of the weekend reading through this thing (I'm about halfway through) I can honestly say that this book is well worth the price of admission.

I seem to do these things on a long enough time schedule to forget absolutely everything, but frequently enough that it's still a pain to look them up and chase down the SQL esoterica needed to accomplish a task. So, here are some canned recipes for manipulating things in PostGIS's flavor of SQL.

  1. Adding an SRID to an existing geometry column in a database. (Something about getting worm points out of my code to the various database intermediaries appears to drop the SRID in the process...)

    1. SELECT UpdateGeometrySRID('Schema Name', 'mytable', 'the_geom', newSRID) ;

      (from https://gis.stackexchange.com/questions/34612/how-to-change-the-srid-of-exisisting-data-in-postgis)

  2. Clip a piece of vector geometry with a polygon and write the result to a new table. (Useful to clip worm point or other geometry to the buffered outline of our tri-state area. The SRID of the worm points is being converted on the fly in the WHERE clause.)

    1. CREATE TABLE clipped_app_basin_bga_worm_points AS
SELECT tp.*
FROM public.test_points AS tp, public.buffered_dissolved_ny_pa_wv AS poly
WHERE ST_Intersects(ST_Transform([tp.pt](http://tp.pt),4326),poly.geom);

      (modified from section 4.7.2.5 of http://postgis.net/docs/using_postgis_dbmanagement.html#examples_spatial_sql)

  3. Add a geometry column to a table.

    1. SELECT AddGeometryColumn ( 'MergedBouguerAnomalies','geom', 4326, 'POINT', 2); Or, if that causes trouble with LinestringZM:
    2. ALTER TABLE some_table ADD COLUMN geom geometry(Point,4326);

  4. And then populate the geometry from latitude and longitude fields. N.B. the double quotes needed for PostGIS SQL dialect...

    1. UPDATE adk_merged_eqs SET geom = ST_SetSRID(ST_MakePoint("_Longitude_","_Latitude_"),4326);

  5. Copy geometry from one table to another (useful in merging gravity measurements from different sources). This is also an example of extracting a POINT from a MULTIPOINT via ST_Dump.

    1. INSERT INTO merged_bga (geom, bga)
      SELECT (ST_Dump("PACES_buffered_grav".geom)).geom AS the_POINT_geom, "PACES_buffered_grav".cbanom267 AS the_BGA
      FROM
        public."PACES_buffered_grav";

  6. Clip vector (multi)linestring geometry with a polygon and write the result to a new table. (Useful to clip worm lines to the buffered outline of our tri-state area. The SRID of the worm points is being converted on the fly in the ST_Intersection call.)

    1. CREATE TABLE clipped_app_basin_bga_worms AS

      SELECT ST_Intersection(ST_Transform(tl.geom,4326),poly.geom)

      FROM public.test AS tl, public.buffered_dissolved_ny_pa_wv AS poly;

      (modified from the example in http://postgis.net/docs/ST_Intersection.html)

    2. CREATE TABLE clipped_app_basin_merged_eqs AS

      SELECT tl.*

      FROM merged_eqs AS tl, public.buffered_dissolved_ny_pa_wv AS poly

      WHERE ST_Intersects(ST_Transform(tl.wkb_geometry,4326),poly.geom) (this one preserves all of the fields of the database, not just the geometry...)

  7. A geometry nearby query:

    1. CREATE TABLE worm_points_near_eqs_5km AS SELECT DISTINCT wp.* FROM clipped_app_basin_bga_worm_points AS wp, clipped_app_basin_eqs_ta_to_oct2014 AS eqs WHERE ST_DWithin([wp.pt](http://wp.pt),ST_Transform(eqs.wkb_geometry,32618),5000.0)

  8. Patch up a table (created by a query?) that is missing a Primary Key:

    1. ALTER TABLE worm_points_near_eqs_5km ADD PRIMARY KEY (id);

  9. Build the all-worms-near-earthquakes table. (The conceptual problem I was having was that I didn't understand that joining a table with a filtered version of itself might be useful...)

    CREATE TABLE "ravatAppBasinPSG1250_worms_near_eqs_5km" AS
    SELECT tlp.*
    FROM "ravatAppBasinPSG1250_levels_points" AS tlp
    JOIN
    (SELECT DISTINCT ws.worm_level_id AS wlid, ws.worm_seg_id AS wsi
    FROM "ravatAppBasinPSG1250_levels_points" AS ws, clipped_ta_neic_eqs_app_basin_through_may_2015 AS eqs
    WHERE ST_DWithin(ws.line_segmt,ST_Transform(eqs.geom,32618),5000.0) AND (eqs.bix_potential_blasts=FALSE)) AS filtered_worm_segs
    ON  (wlid = tlp.worm_level_id) AND (wsi = tlp.worm_seg_id)

  10. Construct a polygon for the "slim bounding box" of a cross-section line: 1. 

    UPDATE adk_section_table SET geom = ST_MakePolygon(ST_GeomFromText('LINESTRING(248896 5126903, 903829 4895705, 903829 4895704, 248896 5126902, 248896 5126903)'));
2.  INSERT INTO ont_section_table (geom)  VALUES (ST_MakePolygon(ST_GeomFromEWKT('SRID=32618;LINESTRING(170108 4811987, 402925 4813555, 402925 4803555, 170108 4801987, 170108 4811987)')));
  1. Build a 3D spatial index on a field 1. CREATE INDEX idx_adk_2500_points_pt ON "ADK_2500_points" USING gist(pt gist_geometry_ops_nd);
  2. Build a slim bounding box intersection for constructing a cross-section showing worm points:
1.  <pre>CREATE TABLE adk_b_b_prime_1000_psg_intersections AS</pre>

    <pre>SELECT DISTINCT p.geom AS pt, p.x AS x, p.y AS y, p.z AS z, p.height AS layer, p.grad AS grad, p.worm_point_id AS id</pre>

    <pre>FROM "test_psg_points" AS p , adk_section_table AS poly</pre>

    <pre>WHERE ST_Intersects(p.geom,poly.geom);</pre>

    <pre>ALTER TABLE adk_b_b_prime_1000_psg_intersections ADD PRIMARY KEY (id);</pre>

    <pre>SELECT UpdateGeometrySRID('public', 'adk_b_b_prime_1000_psg_intersections', 'pt', 32618) ;</pre>

  1. Build a table of new rows from one table that do not exist in another table. (Bloody tricky. SQL. Grrrr.) 1. CREATE TABLE new_clipped_ta_events AS

    SELECT t.*

    FROM

       "clipped_app_basin_merged_eqs" AS m

    RIGHT OUTER JOIN

       "clipped_ta_through_march_2015" as t

    ON

       (m._catalog_ = t."_Catalog_") AND (m._eventid_ = t."_EventID_")

    WHERE m.eventid IS Null

2.  <pre>(Then insert the result):</pre>

    INSERT INTO clipped_app_basin_merged_eqs

    (wkb_geometry,

    _eventid_,

    _time_,

    _latitude_,

    _longitude_,

    _depth_km_,

    _author_,

    _catalog_,

    _contributor_,

    _contributorid_,

    _magtype_,

    _magnitude_,

    _magauthor_,

    _eventlocationname

    )

    SELECT

    ST_SetSRID(ST_MakePoint(ST_X(n.geom),ST_Y(n.geom),n."_Depth_km_"),4326),

    n."_EventID_",

    n."_Time_",

    n."_Latitude_",

    n."_Longitude_",

    n."_Depth_km_",

    n."_Author_",

    n."_Catalog_",

    n."_Contributor_",

    n."_ContributorID_",

    n."_MagType_",

    n."_Magnitude_",

    n."_MagAuthor_",

    n."_EventLocationName"

    FROM new_clipped_ta_events AS n

  1. Worms to Max Grad: 1. CREATE TABLE "ADK_2500_levels_points_to_max_grad" AS

    SELECT DISTINCT lower.*

    FROM

    (

      SELECT a.worm_level_id,a.line_segmt,a.line_grad

  FROM "ADK_2500_levels_points" AS a

    ) as upper

    JOIN

    (

      SELECT b.*

  FROM "ADK_2500_levels_points" AS b

    ) AS lower

    ON ( (lower.worm_level_id = upper.worm_level_id+1) AND (lower.line_grad > upper.line_grad) )

    WHERE ST_3DDWithin(upper.line_segmt,lower.line_segmt,5000.0)

  2. In preparing final Phase1 data for upload to the NGDS/GDR, it became clear that having all tables under the PUBLIC schema is a mistake. The following incantation fixes that: 1. CREATE SCHEMA gpfa_ab_final_earthquake_data

    ALTER TABLE public.clipped_ta_neic_eqs_app_basin_through_may_2015 SET SCHEMA gpfa_ab_final_earthquake_data

    16. To really speed up spatial queries, it's not enough to merely index the spatial components. You also need to "cluster" them. Here's an example:

            a. CREATE INDEX "ADKFinalBGA2500_to_max_grad_levels_points_idx" ON "ADKFinalBGA2500_to_max_grad_levels_points" USING GIST (line_segmt);

            CLUSTER "ADKFinalBGA2500_to_max_grad_levels_points" USING "ADKFinalBGA2500_to_max_grad_levels_points_idx";

        b. You then need to "VACUUM ANALYZE" the table.

        c. For the ADK mag worm earthquake proximity  to extended worms query, this trick took the runtimes from about 20 minutes to 3.6 seconds~~!~~

This is not really PostGIS per se, but here is the incantation to bring a spatialite database table over to PostGIS:

  ogr2ogr -f "PostgreSQL" PG:"dbname=frank" ANFEventsThrough10_2014.sqlite  ta_anf_events_oct2014

Another incantation to bring a csv file into a Postgres database using ogr2ogr.

ogr2ogr -f "PostgreSQL" PG:"port=5433 user=frank password=mypasshere dbname=frank" -nln "adk_bga_euler_new" ADK_BGA_Euler_New.csv

WARNING! I can't figure out how to have the field type autodetected from that command line invocation. So, I resorted to a foo.csvt file matching with the foo.csv file with contents of the form:

"Real","Real","Real","Real","Real","Real","Real","Real","Real","Real","Real","Integer"

Very usefully, a method exists to hook up a Postgres/PostGIS database as a virtual thingy in sqlite/spatialite. This is the only way I could figure out how to reliably get (complicated structurally) worm data from PostGIS to spatialite: https://www.gaia-gis.it/fossil/virtualpg/wiki?name=tutorial This works.

There is an attached file on this wiki page that contains the full blown spatialite/sql script to migrate a levels_points linking table from PostGIS to spatialite. A different attached file does the whole job for all 3 interlinked tables, not just the levels_points table like in the first example. The text follows...

/* spatialite sql code that brings a PostGIS *_levels_points table over to spatialite, and keeps the critical structures intact. */

/* Now brings the levels and points tables too */

/* Load the PostGIS access extension */

SELECT load_extension('mod_virtualpg');

/* Stop it whingeing about violating foreign key constraints as we build it */

PRAGMA foreign_keys=OFF;

/* Create levels table with the right structures */

CREATE TABLE EarBga5km_levels(

worm_level_id INT,

level REAL,

PRIMARY KEY(worm_level_id)

);

/* Connect to the PostGIS version */

CREATE VIRTUAL TABLE pg_bga_levels USING VirtualPostgres ('host=localhost port=5432 dbname=frank user=frank', public, EarBga5km_levels);

/* Populate the rows of the local table with the PostGIS data */

INSERT INTO EarBga5km_levels SELECT * from pg_bga_levels;

/* Close off the PostGIS connection to avoid potential problems */

DROP TABLE pg_bga_levels;

/* That's it for the levels table. Now to the points table... */

CREATE TABLE EarBga5km_points(

worm_point_id INT,

vtk_id INT,

x REAL,

y REAL,

z REAL,

grad REAL,

height REAL,

pt NUM,

wgs84_pt NUM,

PRIMARY KEY(worm_point_id)

);

/* Connect to the PostGIS version */

CREATE VIRTUAL TABLE pg_bga_points USING VirtualPostgres ('host=localhost port=5432 dbname=frank user=frank', public, EarBga5km_points);

/* Populate the rows of the local table with the PostGIS data */

INSERT INTO EarBga5km_points SELECT * from pg_bga_points;

/* Close off the PostGIS connection to avoid potential problems */

DROP TABLE pg_bga_points;

/* Now, munge the geometry fields into something spatialite understands */

UPDATE EarBga5km_points SET pt = GeomFromEWKB(pt);

SELECT RecoverGeometryColumn( 'EarBga5km_points', 'pt', 32737, 'POINT', 'XYZM');

UPDATE EarBga5km_points SET wgs84_pt = GeomFromEWKB(wgs84_pt);

SELECT RecoverGeometryColumn( 'EarBga5km_points', 'wgs84_pt', 4326, 'POINT', 'XYZM');

/* That's it for the points table. Now to the levels_points table... */

/* Create levels_points table with the right structures */

CREATE TABLE EarBga5km_levels_points(

worm_seg_id INT,

seg_sequence_num INT,

line_segmt NUM,

line_grad REAL,

azimuth REAL,

start_point_id INT,

wgs84_line_segmt NUM,

worm_level_id INT,

point_id INT,

PRIMARY KEY(worm_seg_id,worm_level_id,point_id),

FOREIGN KEY(worm_level_id) REFERENCES "EarBga5km_levels_points" (worm_level_id),

FOREIGN KEY(point_id) REFERENCES "EarBga5km_levels_points" (worm_point_id)

);

/* Connect to the PostGIS version */

CREATE VIRTUAL TABLE pg_bga_levels_points USING VirtualPostgres ('host=localhost port=5432 dbname=frank user=frank', public, EarBga5km_levels_points);

/* Stop it whingeing about violating foreign key constraints as we build it */

/* PRAGMA foreign_keys=OFF; */

/* Populate the rows of the local table with the PostGIS data */

INSERT INTO EarBga5km_levels_points SELECT * from pg_bga_levels_points;

/* Close off the PostGIS connection to avoid potential problems */

DROP TABLE pg_bga_levels_points;

/* Now, munge the geometry fields into something spatialite understands */

UPDATE EarBga5km_levels_points SET line_segmt = GeomFromEWKB(line_segmt);

SELECT RecoverGeometryColumn( 'EarBga5km_levels_points', 'line_segmt', 32737, 'LINESTRING', 'XYZM');

UPDATE EarBga5km_levels_points SET wgs84_line_segmt = GeomFromEWKB(wgs84_line_segmt);

SELECT RecoverGeometryColumn( 'EarBga5km_levels_points', 'wgs84_line_segmt', 4326, 'LINESTRING', 'XYZM');

/* And we should have a good levels_points table in spatialite now... */

  1. A proximity query to set the "surveyed" field in a "planned" point database if a point is within 100m of a "as surveyed" table.

UPDATE cornell_grav."CornellDeformedGravStationsWithOwner" AS gp SET surveyed = TRUE FROM cornell_grav."GravSurvey_2018-07-09" AS sp WHERE ST_DWithin(ST_Transform(sp.geom,32618), gp.geom, 100.0)

  1. Add a GeoTIFF raster to a PostGIS database to enable spatial queries of the values (the -t option tiles for performance):

    raster2pgsql -t 200x200 -I -C -s 26918 ./tc_dem_clipped.tif cornell_grav.tc_dem_clipped | psql -d frank

  2. Read raster DEM value (in feet, but converted to meters) at locations of survey points: SELECT id, (ST_Value(rt.rast, 1, pt.geom) * 0.0254 * 12) as RasterValue FROM cornell_grav.tc_dem_clipped_wgs84 AS rt, cornell_grav."GravSurvey_2018-07-16-18-38-38" AS pt WHERE ST_Intersects(pt.geom, rt.rast)

  3. Find duplicate rows in a table, fast: select * from ( SELECT id, ROW_NUMBER() OVER(PARTITION BY column1, column2 ORDER BY id asc) AS Row FROM tbl ) dups where dups.Row > 1

  4. Set a timestamp field from the string in the record: UPDATE cornell_grav."GravSurvey_2018-07-20-17-02-01" AS surv SET timestamp = to_timestamp(surv.created,'yyyy-mm-dd-hh24-mi-ss EDT')

  5. And here's the equivalent for updating the GravityMeasurements table from data off the CG5 meter:

    UPDATE cornell_grav."CG5Measurements" AS cg5 SET cg5_timestamp = to_timestamp(cg5.date ||' '|| cg5.time,'yyyy/mm/dd hh24:mi:ss EDT')

  6. Here's a good place to find out that some of my idioms above are not good ideas: https://wiki.postgresql.org/wiki/Don%27t_Do_This.

Clone this wiki locally