PostGIS - Tipps und Tricks: Unterschied zwischen den Versionen

Aus Geoinformation HSR
Wechseln zu: Navigation, Suche
(PostGIS-Daten laden)
(PostGIS optimieren)
 
(106 dazwischenliegende Versionen von 4 Benutzern werden nicht angezeigt)
Zeile 1: Zeile 1:
Siehe auch:  
+
See also:  
* [[PostGIS]], [[PostGIS-Beispiele]] und [[PostGIS Snippets]]
+
* [[PostGIS Terminal Examples]], [[PostGIS-Beispiele]] und [[PostGIS Snippets]] ([[PostGIS]])
* [[GISpunkt-Seminar_PostGIS]]  
+
* [[PostgreSQL - Tipps und Tricks]] ([[PostgreSQL]])
* [[PostgreSQL]], [[PostgreSQL - Tipps und Tricks]]
+
* Weitere Tipps und Tricks findet man im [http://spatialdbadvisor.com/postgis_tips_tricks SpatialDBAdvisor Blog für PostGIS]
  
  
== Tutorial zum Erstellen einer PostGIS-Datenbank ==
+
== PostGIS-Datenbank mit Geometrie-Tabellen erzeugen ==
  
 
Quelle: [http://www.bostongis.com/?content_name=postgis_tut01 BostonGIS].
 
Quelle: [http://www.bostongis.com/?content_name=postgis_tut01 BostonGIS].
Zeile 11: Zeile 11:
 
Erläuterungen: "shell>" ist der Prompt einer DOS- oder Linux-Shell. "gisdb=#" ist der Prompt von psql, eingeloggt als gisdb-User.
 
Erläuterungen: "shell>" ist der Prompt einer DOS- oder Linux-Shell. "gisdb=#" ist der Prompt von psql, eingeloggt als gisdb-User.
  
;PostGIS-Versionen kontrollieren:
+
;PostgreSQL und PostGIS-Versionen kontrollieren:
 
* PostgreSQL-Version: SELECT version();
 
* PostgreSQL-Version: SELECT version();
 
* PostGIS-Version:  
 
* PostGIS-Version:  
Zeile 17: Zeile 17:
 
* sonst template kontrollieren?
 
* sonst template kontrollieren?
  
;PostGIS-Datenbank erzeugen (Name gisdb):
+
;PostGIS-Extension installieren
* Zur Vorbereitung: User 'gisdb' erzeugen:
+
   gisb=# CREATE EXTENSION postgis;
  shell> createuser -a -D gisdb
 
    oder
 
   shell> createuser -s -D -U postgres gisdb
 
  
* Datenbank 'gisdb' mit Template 'template_postgis' und User 'gisdb' erzeugen:
+
;PostGIS-Tabelle erstellen (mit GEOMETRY-Typ und "type modifier (typmod))):
  shell> createdb -O gisdb -T template_postgis -U postgres gisdb description
+
   CREATE TABLE geofoo (
 
+
    id int4 PRIMARY KEY,
* Session und Grants festlegen:
+
    name text,
  shell> psql -U postgres template1
+
    geom geometry('POINT', 21781)
  template1=# \c gisdb
+
   )
  You are now connected to database "gisdb".
 
  gisdb=# GRANT ALL ON TABLE geometry_columns TO postgres,public;
 
  GRANT
 
  gisdb=# GRANT ALL ON TABLE spatial_ref_sys TO postgres,public;
 
  GRANT
 
 
 
* DB-Projektion (default SRID=-1) vormerken:
 
  gisdb=# select srid, srtext, proj4text from spatial_ref_sys where srtext ILIKE '%Massachusetts%';
 
 
 
;Download vornehmen:
 
* Downloaden: ftp://data.massgis.state.ma.us/pub/shape/state/towns.exe
 
* Auspacken und überprüfen:
 
  shell> unzip towns.exe
 
  shell> ogrinfo -so -al towns_poly.shp
 
 
 
;Shapefile in PostGIS laden:
 
Shapefile nach SQL konvertieren: (-S erzeugt 'single' Polygone (was shapefiles sind) und keine Multipolygone):
 
   shell> shp2pgsql -s 26986 -S TOWNS_POLY towns >towns_poly.sql
 
  oder
 
  shell> shp2pgsql -s 4326 -W latin1 Kantone_WGS84 kantone >kantone.sql
 
 
 
;SQL-Daten in PostGIS laden: ('2>towns_psql_error.txt' leitet Error-Meldungen in eine Datei um):
 
   shell> psql -d gisdb -h localhost -U gisdb -f towns_poly.sql >towns_psql_log.txt 2>towns_psql_error.txt
 
  oder
 
  shell> psql -U postgres -d gisdb2 -f kantone.sql
 
 
 
Beim Import können u.a. Topologie-Fehler zum Vorschein kommen ("Error: geometry contains non-closed rings"). Solche Geometrien werden nicht importiert.
 
  
 
;PostGIS-Index erzeugen und Installation überprüfen:  
 
;PostGIS-Index erzeugen und Installation überprüfen:  
Zeile 62: Zeile 32:
 
   gisdb=# CREATE INDEX idx_towns_the_geom ON towns USING gist(the_geom);
 
   gisdb=# CREATE INDEX idx_towns_the_geom ON towns USING gist(the_geom);
 
   gisdb=# CREATE INDEX idx_towns_town ON towns USING btree(town);
 
   gisdb=# CREATE INDEX idx_towns_town ON towns USING btree(town);
 
PostGIS-DB-Installation überprüfen:
 
  gisdb=# SELECT extent(the_geom) FROM towns WHERE town = 'BOSTON';
 
  gisdb=# SELECT area(geomunion(the_geom)) FROM towns WHERE town = 'BOSTON';
 
  gisdb=# SELECT town,astext(transform(the_geom,4326)) FROM towns WHERE town = 'BOSTON';
 
  gisdb=# SELECT count(town) FROM towns WHERE town = 'BOSTON';
 
  gisdb=# SELECT st_askml(geomunion(transform(the_geom,4326))) AS the_geom FROM towns WHERE town = 'BOSTON';
 
 
;PostGIS-DB/-Tabelle löschen:
 
* PostGIS-DB:
 
  shell> dropdb gisdb
 
* Tabelle löschen inkl. aller Einträge in den "Systemtabellen":
 
  abhängigkeiten in den 'Ref-System'-Tabellen?
 
  
 
== Shapefiles in PostGIS importieren (shp2pgsql) ==
 
== Shapefiles in PostGIS importieren (shp2pgsql) ==
Zeile 89: Zeile 46:
  
 
* Mittels (mit PostGIS ausgeliefertem) Kommandozeilen-Tool 'shp2pgsql'.  
 
* Mittels (mit PostGIS ausgeliefertem) Kommandozeilen-Tool 'shp2pgsql'.  
* Mittels INSERT und WKT oder Konstruktoren.
+
* Mittels INSERT und WKT oder Konstruktoren (vgl. Snippet "Creating Geometry (...) Types" unten).
  
Allowed and non-functioning constructors for geometry types (Varianten für geometry types):
+
=== Creating Geometry and Geography Types ===
 +
Allowed and non-functioning constructors for creating '''geometry''' type POINT (lon lat, +/-180 +/-90):
 
   /*0*/ SELECT ST_AsEWKT(ST_GeomFromText('POINT(-71.06 42.28)'))          -- '''Preferred simplest text form without SRID'''
 
   /*0*/ SELECT ST_AsEWKT(ST_GeomFromText('POINT(-71.06 42.28)'))          -- '''Preferred simplest text form without SRID'''
 
   /*1*/ SELECT ST_AsEWKT(ST_GeomFromText('POINT(-71.06 42.28)', 4326))    -- '''Preferred for text form with SRID'''
 
   /*1*/ SELECT ST_AsEWKT(ST_GeomFromText('POINT(-71.06 42.28)', 4326))    -- '''Preferred for text form with SRID'''
Zeile 98: Zeile 56:
 
   /*4*/ SELECT ST_AsEWKT('SRID=4326;POINT(-71.06 42.28)'::geometry)        -- With cast; prefer *1*/ST_GeomFromText
 
   /*4*/ SELECT ST_AsEWKT('SRID=4326;POINT(-71.06 42.28)'::geometry)        -- With cast; prefer *1*/ST_GeomFromText
 
   /*5*/ SELECT ST_AsEWKT(ST_SetSRID('POINT(-71.06 42.28)'::geometry,4326)) -- With cast; prefer *1*/ST_GeomFromText
 
   /*5*/ SELECT ST_AsEWKT(ST_SetSRID('POINT(-71.06 42.28)'::geometry,4326)) -- With cast; prefer *1*/ST_GeomFromText
   /*6*/ SELECT ST_AsEWKT(ST_MakePoint(-71.06, 42.28, 4326))               -- '''Preferred symbolic form''' (Hint: returns WKT, not EWKT)
+
   /*6*/ SELECT ST_AsEWKT(ST_SetSRID(ST_MakePoint(-71.06, 42.28),4326))     -- '''Preferred symbolic form with EWKT'''
  
Allowed and non-functioning constructors for geography types (Varianten für geography types):
+
Allowed and non-functioning constructors for creating '''geography''' types (Varianten für geography types):
 
   /*1*/  -- SELECT ST_AsEWKT(ST_GeogFromText('POINT(-71.06 42.28)', 4326))      -- ERROR: unknown ST_GeogFromText()
 
   /*1*/  -- SELECT ST_AsEWKT(ST_GeogFromText('POINT(-71.06 42.28)', 4326))      -- ERROR: unknown ST_GeogFromText()
 
   /*2a*/ -- SELECT ST_AsEWKT(ST_GeogFromText('SRID=4326;POINT(-71.06 42.28)'))  -- ERROR: no ST_AsEWKT: why not?
 
   /*2a*/ -- SELECT ST_AsEWKT(ST_GeogFromText('SRID=4326;POINT(-71.06 42.28)'))  -- ERROR: no ST_AsEWKT: why not?
   /*2b*/ SELECT ST_AsText(ST_GeogFromText('SRID=4326;POINT(-71.06 42.28)'))      -- Preferred with srid 'inline'; (Hint: ST_AsEWKT doesn't work)
+
   /*2b*/ SELECT ST_AsText(ST_GeogFromText('SRID=4326;POINT(-71.06 42.28)'))      -- '''Preferred with srid 'inline';''' (Hint: ST_AsEWKT doesn't work)
 
   /*3*/  -- SELECT ST_AsText(ST_GeogFromEWKT('SRID=4326;POINT(-71.06 42.28)'))  -- ERROR: no ST_GeogFromEWKT()
 
   /*3*/  -- SELECT ST_AsText(ST_GeogFromEWKT('SRID=4326;POINT(-71.06 42.28)'))  -- ERROR: no ST_GeogFromEWKT()
 
   /*4*/  SELECT ST_AsText('SRID=4326;POINT(-71.06 42.28)'::geography)            -- With cast; prefer *1*; (Hint: ST_AsEWKT doesn't work)
 
   /*4*/  SELECT ST_AsText('SRID=4326;POINT(-71.06 42.28)'::geography)            -- With cast; prefer *1*; (Hint: ST_AsEWKT doesn't work)
 
   /*5*/  --- SELECT ST_AsText(ST_SetSRID('POINT(-71.06 42.28)'::geography,4326)) -- ERROR: ok SetSRID w. geography unnecessary
 
   /*5*/  --- SELECT ST_AsText(ST_SetSRID('POINT(-71.06 42.28)'::geography,4326)) -- ERROR: ok SetSRID w. geography unnecessary
   /*6*/  -- Equivalent to ST_MakePoint() for geography types missing.            -- .
+
   /*6*/  SELECT ST_MakePoint(-71.06, 42.28)                                     -- for geography types missing
 +
 
 +
=== Create Polygon given Bounding Box (BBox) Coordinates ===
 +
 
 +
Try this in http://labs.geometa.info/postgisterminal :
 +
 
 +
  SELECT ST_AsText(ST_Transform(ST_MakeEnvelope(8.795611, 46.872886, 9.674135, 47.675419, 4326), 3857)) AS geom, '#' AS label
 +
  SELECT ST_AsText(ST_Transform(ST_SetSRID(ST_Envelope('LINESTRING(8.795611 46.872886, 9.674135 47.675419)'::geometry),4326), 3857)) AS geom, '#' AS label
 +
  SELECT ST_AsText(ST_Transform(ST_SetSRID('BOX(8.795611 46.872886, 9.674135 47.675419)'::box2d, 4326), 3857)) AS geom, '#' AS label
 +
  SELECT ST_AsText(ST_Transform(ST_SetSRID('BOX3D(8.795611 46.872886, 9.674135 47.675419)'::box3d, 4326), 3857)) AS geom, '#' AS label
 +
  SELECT ST_AsText(ST_SetSRID(ST_Extent(way), 4326)) AS geom FROM osm_point
 +
 
 +
== Curved Geometries ==
 +
 
 +
Curved geometries in PostGIS with notes for QGIS (and OGR).
 +
 
 +
Introduction: See '''GeoPackage Spec. (http://www.geopackage.org/spec110/#sfsql_intro)''' und '''PostGIS 2.5 Manual (https://postgis.net/docs/manual-2.5/'''. Citing PostGIS 2.5 Manual "4.1.3 SQL-MM Part 3": The SQL Multimedia Applications Spatial specification extends the simple features for SQL spec by defining a number of circularly interpolated curves. The SQL-MM definitions include 3DM, 3DZ and 4D coordinates, but do not allow the embedding of SRID information. The Well-Known Text extensions are not yet fully supported.". This diagram is interestingly showing the SQL-MM geometry types: http://www.dpi.inpe.br/terralib5/wiki/doku.php?id=wiki:documentation:devguide:geometry_module :
 +
 
 +
[[Datei:Iso_sql_mm_geometry_type_hierarchy.png]]
 +
 
 +
These are the curved types (from the PostGIS 2.5 Manual):
 +
* CIRCULARSTRING is the basic curve type, similar to a LINESTRING. A single segment has second point as point on the arc (there are also closed circles, chained arcs)
 +
* COMPOUNDCURVE is a single, continuous curve that has both curved (circular) segments and linear segments.
 +
* CURVEPOLYGON is just like a polygon. The difference is that a outer/inner ring can take the form of a circular string, linear string or compound string.
 +
* MULTICURVE is a collection of curves, which can include linear strings, circular strings or compound strings.
 +
* MULTISURFACE is a collection of surfaces, which can be (linear) polygons or curve polygons.
 +
 
 +
Note that MultiCurve, MultiSurface are called abstract (and reserved) in the OGC model. Astonishingly in PostGIS, GeoPackage and QGIS it's used anyway.
 +
 
 +
QGIS: Supports these (CircularString, CompoundCurve, CurvePolygon, MultiCurve, MultiSurface). See https://www.qgis.ch/de/ressourcen/anwendertreffen/2015/curved-geometries-interlis-with-ogr-and-qgis .
 +
 
 +
OGR: Some functions to convert curves/arcs to linstrings.
 +
 
 +
Examples: (Note one could omit LINESTRING in 'LINESTRING(1 0, 0 1))' and also POLYGON here):
 +
    SELECT ST_AsText('CIRCULARSTRING(0 0, 1 1, 1 0)'::geometry);
 +
    SELECT ST_AsText('COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 1 0),LINESTRING(1 0, 0 1))'::geometry);
 +
    SELECT ST_AsText('CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0,2 0, 2 1,2 3,4 3),LINESTRING(4 3, 4 5, 1 4, 0 0)), CIRCULARSTRING(1.7 1, 1.4 0.4, 1.6 0.4, 1.6 0.5, 1.7 1))'::geometry);
 +
    SELECT ST_AsText('MULTICURVE((0 0, 5 5),CIRCULARSTRING(4 0, 4 4, 8 4))'::geometry);
 +
    SELECT ST_AsText('MULTISURFACE(CURVEPOLYGON(CIRCULARSTRING(0 0, 4 0, 4 4, 0 4, 0 0),LINESTRING(1 1, 3 3, 3 1, 1 1)),POLYGON((10 10, 14 12, 11 10, 10 10),(11 11, 11.5 11, 11 11.5, 11 11)))'::geometry);
 +
 
 +
 
 +
Functions:
 +
* ST_IsCollection returns TRUE if the geometry type of the argument is either a GEOMETRYCOLLECTION, MULTI{POINT,POLYGON,LINESTRING,CURVE,SURFACE}, or a COMPOUNDCURVE. NOTE: This PostGIS-Doc. should also mention CURVEPOLYGON!
 +
* ST_HasArc: e.g. show if it's an art aber breaking up a compound curve into its segments
 +
    SELECT ST_AsEWKT(a.geom), ST_HasArc(a.geom)
 +
    FROM (SELECT (ST_Dump(p_geom)).geom AS geom
 +
      FROM (SELECT ST_GeomFromEWKT('COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 1 0),(1 0, 0 1))') AS p_geom) AS b
 +
    ) AS a;
 +
    -- "CIRCULARSTRING(0 0,1 1,1 0)"    true
 +
    -- "LINESTRING(1 0,0 1)"           false
 +
* ST_ForceCurve — Upcast a geometry into its curved type, if applicable.
 +
    SELECT ST_AsText(
 +
        ST_ForceCurve(
 +
          'POLYGON((0 0, 5 0, 0 5, 0 0),(1 1, 1 3, 3 1, 1 1))'::geometry
 +
        )
 +
    );
 +
    -- "CURVEPOLYGON((0 0,5 0,0 5,0 0),(1 1,1 3,3 1,1 1))"
 +
* ST_LineToCurve — Converts a LINESTRING/POLYGON to a CIRCULARSTRING, CURVEPOLYGON. If is not curved enough to clearly represent a curve the function will return the same input geometry.
 +
* ST_CurveToLine(geometry curveGeom, float tolerance, integer tolerance_type, integer flags) — Converts a CIRCULARSTRING/CURVEPOLYGON/MULTISURFACE to a LINESTRING/POLYGON/MULTIPOLYGON
 +
 
 +
 
 +
== PostGIS-Daten exportieren ==
 +
 
 +
=== Export a table with a geometry attribute to GeoJSON ===
 +
 
 +
Sources:
 +
* http://blog.cleverelephant.ca/2019/03/geojson.html
 +
* https://gis.stackexchange.com/questions/112057/sql-query-to-have-a-complete-geojson-feature-from-postgis (validated 2018-10-28 SK)
 +
* (Other, older source: http://www.postgresonline.com/journal/archives/267-Creating-GeoJSON-Feature-Collections-with-JSON-and-PostGIS-functions.html )
 +
 
 +
Query:
 +
  :: File rowjsonb_to_geojson.sql ::
 +
  /* ------------------------------------------------------------------
 +
  Usage :
 +
    SELECT rowjsonb_to_geojson(to_jsonb(mytable.*)) FROM mytable;
 +
  Parameters:
 +
    JSONB - jsonb data type (aka one row)
 +
    geom_column - name of geometry column (default 'geom')
 +
  Credits: P. Ramsey, http://blog.cleverelephant.ca/2019/03/geojson.html
 +
  ------------------------------------------------------------------ */
 +
  CREATE OR REPLACE FUNCTION rowjsonb_to_geojson(
 +
    rowjsonb JSONB,
 +
    geom_column TEXT DEFAULT 'geom')
 +
  RETURNS TEXT AS $$
 +
  DECLARE
 +
    json_props jsonb;
 +
    json_geom jsonb;
 +
    json_type jsonb;
 +
  BEGIN
 +
    IF NOT rowjsonb ? geom_column THEN
 +
      RAISE EXCEPTION 'geometry column ''%'' is missing', geom_column;
 +
    END IF;
 +
    json_geom := ST_AsGeoJSON((rowjsonb ->> geom_column)::geometry)::jsonb;
 +
    json_geom := jsonb_build_object('geometry', json_geom);
 +
    json_props := jsonb_build_object('properties', rowjsonb - geom_column);
 +
    json_type := jsonb_build_object('type', 'Feature');
 +
    return (json_type || json_geom || json_props)::text;
 +
  END;
 +
  $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT;
  
== PostGIS-Daten darstellen ==
+
(Query - older solution):
 +
  :: File togeojson.sql ::
 +
  /* ---------------------------
 +
  Convert following query output to GeoJSON:
 +
    SELECT
 +
      st_transform(geom,4326) AS geom, -- geometry
 +
      gid as id, "name", "level", pop_class -- "properties"
 +
    FROM orte
 +
    LIMIT 10
 +
  --------------------------- */
 +
  SELECT jsonb_build_object(
 +
    'type',    'FeatureCollection',
 +
    'features', jsonb_agg(feature)
 +
  ) AS geojson
 +
  FROM (
 +
    SELECT jsonb_build_object(
 +
      'type',      'Feature',
 +
      --'id',        id, -- 'id' is ignored e.g. by geojson.io
 +
      'geometry',  st_asgeojson(geom)::jsonb,
 +
      'properties', to_jsonb(row) - 'geom' -- removing elements with '-'
 +
    ) AS feature
 +
    FROM ( -- input_table (beware columns named like keywords)
 +
      SELECT
 +
        st_transform(geom,4326) AS geom,
 +
        gid as id, "name", "level", pop_class
 +
      FROM orte
 +
      LIMIT 10
 +
    ) AS row
 +
  ) AS features;
  
=== Übersicht ===
+
Call it e.g. using psql like this:
* Lokal, mittels Desktop GIS:
+
  % psql -U <username> -d <dbname> -f togeojson.sql -t -o out.geojson
** [[QGIS]]
+
Ev. beautify output using geojson.io or http://geojsonlint.com/
** [[OpenJUMP]]
 
** GDV Spatial Commander [http://www.gdv.com/down/scommander.php]
 
** gvSIG [http://www.gvsig.gva.es/index.php?id=gvsig&L=2] gvSIG Wiki [http://webmapping.info/mwgvsig/index.php?title=Hauptseite]
 
** etc.
 
* über Webapplikationen ([[WMS]]):
 
** [[GeoServer]]
 
** [[UMN MapServer]]
 
** [[ArcGIS]] Server
 
  
=== PostGIS und KML ===
+
=== Export a table with a geometry attribute to KML ===
  
 
PostGIS kennt die Funktion ST_AsKML(), die eine Geometrie KML-konform ausgibt. Allerdings werden nur diejenigen KML-Elemente der Geometrie selbst ausgegeben. Kopf- und Fusszeile des KML-Dokuments muss man selber hinzufügen. Dazu erzeugen wir nun eine Funktion askmldoc, die ein gültiges KML-Dokument ausgibt:
 
PostGIS kennt die Funktion ST_AsKML(), die eine Geometrie KML-konform ausgibt. Allerdings werden nur diejenigen KML-Elemente der Geometrie selbst ausgegeben. Kopf- und Fusszeile des KML-Dokuments muss man selber hinzufügen. Dazu erzeugen wir nun eine Funktion askmldoc, die ein gültiges KML-Dokument ausgibt:
Zeile 159: Zeile 235:
  
 
  % psql -A -t -d gisdb -c "SELECT askmldoc('MyTable', 'A Comment...', the_geom) FROM mytable WHERE gid=1;" -o mykmlfile.kml
 
  % psql -A -t -d gisdb -c "SELECT askmldoc('MyTable', 'A Comment...', the_geom) FROM mytable WHERE gid=1;" -o mykmlfile.kml
 +
 +
 +
== PostGIS-Daten darstellen ==
 +
 +
Übersicht:
 +
* Lokal, mittels Desktop GIS:
 +
** [[QGIS]]
 +
** [[OpenJUMP]]
 +
** GDV Spatial Commander [http://www.gdv.com/down/scommander.php]
 +
** gvSIG [http://www.gvsig.gva.es/index.php?id=gvsig&L=2] gvSIG Wiki [http://webmapping.info/mwgvsig/index.php?title=Hauptseite]
 +
** etc.
 +
* über Webapplikationen ([[WMS]]):
 +
** [[GeoServer]]
 +
** [[UMN MapServer]]
 +
** [[ArcGIS]] Server
  
 
=== PostGIS und Google Earth ===  
 
=== PostGIS und Google Earth ===  
Zeile 192: Zeile 283:
 
* Beispiel für gisdb:
 
* Beispiel für gisdb:
 
   SELECT ST_AsKML(geomunion(transform(the_geom,4326))) as the_geom from towns where town='BOSTON';
 
   SELECT ST_AsKML(geomunion(transform(the_geom,4326))) as the_geom from towns where town='BOSTON';
 +
 +
== Spatial Relationships and ST_Relate ==
 +
 +
ST_Relate is based on '''Dimensionally Extended Nine-Intersection Model (DE-9IM) (or Clementini-Matrix)'''. It covers all known spatial relationships. The most typical spatial relationsships (and it's opposites) got own functions, like:
 +
* ST_Within != ST_Contains (oder ST_Covers (*))
 +
* ST_Covers != ST_CoveredBy
 +
* ST_Intersects != ST_Disjoint
 +
 +
(*) Note: Prefer ST_Covers over ST_Contains if lines on boundaries count as „inside“ (Source: Martin Davis‘ blog post: http://lin-ear-th-inking.blogspot.ch/2007/06/subtleties-of-ogc-covers-spatial.html )
 +
 +
 +
=== Documentation ===
 +
* Paper [[Media:9dem_springer.pdf|Dimensionally Extended Nine-Intersection Model (DE-9IM)]] - Christian Strobl, DLR, 2007.
 +
* PostGIS documentation ch04: http://postgis.org/documentation/manual-1.5/ch04.html#DE-9IM
 +
* OGC: OpenGIS® Implementation Specification for Geographic information - Simple feature access - Part 1: Common architecture - Extended Nine, 6.1.15.2 The Dimensionally Extended Nine-Intersection Model (DE-9IM)
 +
* Oracle: http://docs.oracle.com/cd/B19306_01/appdev.102/b14255/sdo_operat.htm - Oracle SDO_RELATE() operator ([http://docs.oracle.com/cd/B19306_01/appdev.102/b14255/sdo_intro.htm#i880253])
 +
* MS SQL Server: http://msdn.microsoft.com/en-us/library/bb933915.aspx - SqlServer?2008 STRelate() operator
 +
* Jaspa: http://jaspa.upv.es/jaspa/v0.2.0/manual/html/ST_Relate.html
 +
* Interessant: http://sample.wrd.ir/DE-9IM
 +
 +
=== ST_Relate(geom_a, geom_b, '<<pattern>>') ===
 +
 +
Examples of pattern values:
 +
* ST_Equals:
 +
** 0******** (for Point/Point)
 +
** T*F**FFF* (for any geom?; from Esri Webhelp)
 +
* ST_Contains:
 +
** T*****FF* (for any geom?; from Esri Webhelp)
 +
* ST_Crosses:
 +
** T*T****** (for Point/Line, Point/Area, and Line/Area situations)
 +
** T*****T** (for Line/Point, Area/Point, and Area/Line situations)
 +
** 0******** (for Line/Line situations)
 +
** 0F1FF0102 (for strict Line/Line situation like an 'X')
 +
* ST_Touches (From [http://postgis.org/documentation/manual-1.5/ST_Touches.html]:"The allowable DE-9IM Intersection Matrices for the two geometries are:" ???):
 +
** FT*******  ?
 +
** F**T*****  ?
 +
** F***T****  ?
 +
* ST_Contains(polygon_a, polygon_b) => ST_Relate(polygon_a.geom, table_b.geom, '2********')
 +
* ST_Intersects => 'T********' (or within polygon/polygon: 'T*F**F***' ?)
 +
* ST_Overlaps
 +
** '1********' (for Area/Area)
 +
* ST_Within
 +
** T*F**F*** (for any geom?; from Esri Webhelp)
 +
* disjoint line/point => 'FF*FF****'
 +
* Census blocks and voting district interiors that do not intersect "F" (interior/interior) but have a common linear "1" boundary (boundary/boundary) => (polygonA, polygonB) => 'F***1****'
 +
* What about this topology rule: geometries that neither intersect nor overlap. The only allowed shared points between geometries are boundary-boundary intersections. boundary-interior-boundary (touch interior) intersections are not allowed. (From: Rule.java, Jaspa project [http://jaspa.upv.es/]
 +
 +
=== ST_Relate(geom_a, geom_b) ===
 +
 +
Haben Sie gewusst, dass es auch eine Variante ST_Relate(geomA, geomB) gibt, die die DE-9IM „berechnet“ gegeben einen Geometrie-Input? Hier z.B. zwei sich schneidende Linien (als X) als Input:
 +
 +
  SELECT ST_Relate(PA,PB)
 +
  FROM (
 +
    SELECT ST_GeomFromText('LINESTRING(1 3, 9 6)') AS PA,
 +
          ST_GeomFromText('LINESTRING(3 6, 7 3)') AS PB
 +
    ) AS foo
 +
  => 0F1FF0102
  
 
== PostGIS und Rasterdaten ==
 
== PostGIS und Rasterdaten ==
Zeile 213: Zeile 361:
 
Beim Aufruf von GDAL/OGR wird die GDAL_DATA environment Variable benötigt.
 
Beim Aufruf von GDAL/OGR wird die GDAL_DATA environment Variable benötigt.
  
== Benutzerdefinierte PostgreSQL/PostGIS-Funktionen ==
+
== Erkennen und Eliminieren von Sliver Polygonen und Spikes ==
  
 
Implementiert als sog. 'Stored Procedures'.
 
Implementiert als sog. 'Stored Procedures'.
Zeile 220: Zeile 368:
  
 
=== Eliminate sliver polygons ===
 
=== Eliminate sliver polygons ===
 +
 +
 +
Sources:
 +
* "Spike analyzer and remover"
 +
** PostGIS NormalizeGeometry: https://gasparesganga.com/labs/postgis-normalize-geometry/
 +
** Schmidt, Andreas; Krüger, Nils: spike_analyzer.sql 2009-12-01, Version 1.0, Informatikzentrum Landesverwaltung Baden-Württemberg (IZLBW), url: http://trac.osgeo.org/postgis/wiki/UsersWikiExamplesSpikeAnalyzer.
 +
** Schmidt, Andreas; Krüger, Nils: spikeremover und spikeRemoverCore.sql 2009-10-01, Version 1.0, Informatikzentrum Landesverwaltung Baden-Württemberg (IZLBW), url: http://trac.osgeo.org/postgis/wiki/UsersWikiExamplesSpikeRemover.
 +
* Birgit Laggner und Helge Meyer-Borstel, „Geoprocessing von Massendaten in PostGIS - Probleme und Lösungsansätze“, FOSSGIS 2010, [http://www.fossgis.de/konferenz/2011/programm/track/Vortr%C3%A4ge%20%28GIS%29/170.de.html]
 +
** Komplettverschneidung: http://trac.osgeo.org/postgis/wiki/UsersWikiExamplesOverlayTables2
 +
** Überlappungsbereinigung: http://trac.osgeo.org/postgis/wiki/UsersWikiExamplesPolygonOverlaps
 +
  
 
Given a polygon table that has many small areas and holes. How to remove "small" areas and holes (smaller than a given area in m2)?  
 
Given a polygon table that has many small areas and holes. How to remove "small" areas and holes (smaller than a given area in m2)?  
Zeile 225: Zeile 384:
 
Remarks:  
 
Remarks:  
 
* Similar like the ELIMINATE command in [[ArcGIS]].
 
* Similar like the ELIMINATE command in [[ArcGIS]].
* See also CLEAN topology
+
* See also section [[#Clean topology|Clean topology]] below
  
 
   CREATE OR REPLACE FUNCTION Filter_Rings(geometry,float)
 
   CREATE OR REPLACE FUNCTION Filter_Rings(geometry,float)
Zeile 244: Zeile 403:
 
   $$
 
   $$
 
   LANGUAGE 'sql' IMMUTABLE;
 
   LANGUAGE 'sql' IMMUTABLE;
 +
 +
  CREATE OR REPLACE FUNCTION Filter_Rings(geometry GEOMETRY, param FLOAT)
 +
  RETURNS geometry AS
 +
  $$
 +
  SELECT ST_MakePolygon(c.outer_ring, d.inner_rings) as final_geom
 +
  SELECT * FROM foo WHERE attr > param; -- statt $2 steht hier ein benanntes Fn.-Argument
 +
  ...
 +
  $$
 +
  LANGUAGE 'sql';
  
 
Usage example:  
 
Usage example:  
Zeile 284: Zeile 452:
 
Version:
 
Version:
 
* 1.4 and GEOS 3.1 will bring prepared geometries which will make things faster.
 
* 1.4 and GEOS 3.1 will bring prepared geometries which will make things faster.
 +
 +
Fast map intersection by spatial partioning https://github.com/gbb/fast_map_intersection
  
 
Reorder on disk:
 
Reorder on disk:
Zeile 329: Zeile 499:
 
;Was für Mass-Einheiten (units of measurement) benützen die PostGIS-Funktionen?: Siehe [http://postgis.refractions.net/support/wiki/ PostGIS Wiki].
 
;Was für Mass-Einheiten (units of measurement) benützen die PostGIS-Funktionen?: Siehe [http://postgis.refractions.net/support/wiki/ PostGIS Wiki].
  
;Wie erhält man Distanz-Masse in Meter?: Siehe [http://postgis.refractions.net/support/wiki/ PostGIS Wiki].
+
;Wie erhält man Distanz-Masse in Meter?: Siehe [http://trac.osgeo.org/postgis/wiki/UsersWikiMain PostGIS Wiki].
  
 
;Was sind die Grenzen der 3D-Typen und -Funktionen in PostGIS?: (1) 3D volumetric objects are not supported. (2) 3D non-volume are supported partly -- e.g. a 2D polygon in 3 d space, a line, point in 3d space. (3) Spatial relationships however only consider the spatial component of the x, y plane.
 
;Was sind die Grenzen der 3D-Typen und -Funktionen in PostGIS?: (1) 3D volumetric objects are not supported. (2) 3D non-volume are supported partly -- e.g. a 2D polygon in 3 d space, a line, point in 3d space. (3) Spatial relationships however only consider the spatial component of the x, y plane.
Zeile 337: Zeile 507:
 
== Snippets ==
 
== Snippets ==
  
 +
=== Find centre and radius of a curve defined by three points ===
 +
 +
“COGO: Finding centre and radius of a curve defined by three points (PostGIS)” by Simon Greener:
 +
https://spatialdbadvisor.com/postgis_tips_tricks/278/postgis-cogo-finding-centre-and-radius-of-a-curve-defined-by-three-points
 +
 +
=== Get area of a polygon which has SRID EPSG:4326 (WGS84/GPS) ===
 +
 +
All area functions are squared by definition.
 +
But a sq degree varies in its area depending on latitude.
 +
Transform your data using a SRID referring to an equal-area (or near enough to
 +
equal area for your purposes) projection in the query for area:
 +
 +
  -- Get area from a geometry given the geometry is in 4326.
 +
  % SELECT ST_Area2d(ST_transform(geometry, <SRID>));
  
 
=== Extent of Database ===
 
=== Extent of Database ===
Zeile 395: Zeile 579:
 
   end loop;
 
   end loop;
  
=== Select nearest point on polyline ===
+
=== Find nearest existing point on a linestring given a point (Snapping 1) ===
  
A user clicks somewhere on a polyline and we want to get the closest point that is part of the representation of the line. You can replace 0.01 with some other tolerance number to compensate for the fact that a person probably won't click right on the line.
+
A user clicks somewhere on a polyline and we want to get the closest existing point that is part of the representation of the line. You can replace 0.01 with some other tolerance number to compensate for the fact that a person probably won't click right on the line.
  
 
Ab PostGIS 1.3.5 (wegen generate_series() )?
 
Ab PostGIS 1.3.5 (wegen generate_series() )?
Zeile 406: Zeile 590:
 
   FROM point_table) AS P
 
   FROM point_table) AS P
 
   WHERE ST_Dwithin(P.geom, ST_SetSRID(ST_MakePoint(lon,lat),somesrid), 0.01)
 
   WHERE ST_Dwithin(P.geom, ST_SetSRID(ST_MakePoint(lon,lat),somesrid), 0.01)
 +
 +
=== Find nearest interpolated point on a linestring given a point (Snapping 2) ===
 +
 +
See http://blog.cleverelephant.ca/2008/04/snapping-points-in-postgis.html
 +
 +
=== Find nearest (whole) linestring given a point (Snapping 3) ===
 +
 +
Basic principle given mypos => 'SRID=3005;POINT(1011102 450541)'
 +
 +
Solution using distance within max. radius e.g. in meters (ST_DWithin):
 +
  SELECT geom, name, etc FROM mytable
 +
  WHERE ST_DWithin(mypos, geom, 50) -- max. distance 50m
 +
  ORDER BY ST_Distance(mypos, geom)
 +
  LIMIT 1
 +
 +
and/or solution using 50 nearest objects (kNN index). Note that's returng 50 objects even if they are far away.
 +
  WITH tmp AS (
 +
    SELECT geom, name, ST_Distance(geom, mypos) AS distance
 +
    FROM mytable
 +
    ORDER BY geom <#> mypos -- prefering <#> over <-> when lines are involved.
 +
    LIMIT 50
 +
  )
 +
  SELECT * FROM tmp
 +
  ORDER BY distance
 +
  LIMIT 1
 +
 +
See also the examples at [[PostGIS Terminal Examples]].
 +
 +
 +
=== Find and add nearest hydrant to a parcel ===
 +
 +
Given parcels (or other kind of (point?) feature class) and hydrants (or other kind of (point?) feature class) as input (e.g. imported from external files), find and add closest hydrant to parcel. Output: All parcel attributes incl. geometry, hydrant_id and distance from hydrant to parcel.
 +
 +
This is an example of a lateral join driving a nearest-neighbor distance calculation. Needs PostgreSQL 9.5 / PostGIS 2.x. (Credits: https://twitter.com/pwramsey/status/686980049967443968 )
 +
 +
  SELECT
 +
    -- Take all attributes from parcels incl. pid and geom:
 +
    parcels.*,
 +
    -- Keep the hydrant id around, might be useful later:
 +
    hydrants.id AS hydrant_id,
 +
    -- Calculate distance over the spheroid using geography distance:
 +
    ST_Distance(geography(hydrants.geom), geography(parcels.geom)) AS distance
 +
  FROM
 +
    -- For this data, removing the duplicate parcel geometries and
 +
    -- null attributed parcels cuts the size by 2/3 so a good
 +
    -- performance improvement to filter records here:
 +
    (SELECT DISTINCT ON (geom) *
 +
    FROM parcels_imported
 +
    WHERE pid IS NOT NULL) AS parcels
 +
    -- Cross join takes every combination of join records
 +
    -- without restriction, which is fine since we're limiting
 +
    -- the hydrants side of the join to just one return record
 +
    -- per parcel:
 +
    CROSS JOIN LATERAL
 +
      -- The guts of the query, actually do a KNN query for each
 +
      -- candidate parcel on the other side of the join:
 +
      (SELECT id, geom
 +
      FROM hydrants_imported
 +
      ORDER BY parcels.geom <-> geom -- take here cartesian CRS like webmercator
 +
      LIMIT 1) AS hydrants
 +
    ON TRUE;
  
 
=== PostGIS views ===
 
=== PostGIS views ===
Zeile 419: Zeile 664:
  
 
A3. Grab the key from the underlying geometry table, & use it for a key column in your aggregate view, eg. "create view v1 as select min(key) as key" ...
 
A3. Grab the key from the underlying geometry table, & use it for a key column in your aggregate view, eg. "create view v1 as select min(key) as key" ...
 +
 +
A4. Use window function "SELECT row_number() OVER() As generated_id, ... FROM ..."
  
 
=== Standard direction for polylines ===
 
=== Standard direction for polylines ===
Zeile 585: Zeile 832:
  
 
The split based on GEOMETRYCOLLECTION is that when I first tested the STR indexed ST_Union I found that performance was affected if the union set includes an mpoly.
 
The split based on GEOMETRYCOLLECTION is that when I first tested the STR indexed ST_Union I found that performance was affected if the union set includes an mpoly.
 +
 +
=== Transform coordinates into other reference systems ===
 +
 +
The following statement "abuses" PostGIS to simply transforms the given point coordinates into other reference systems and prints the result to the console. This is useful for example to double-check self coded transformation algorithms.
 +
 +
  SELECT
 +
    ST_AsText(ST_Transform( ST_SetSRID(ST_MakePoint(8.8211, 47.2221), 4326), 4326)) "EPSG:4326",
 +
    ST_AsText(ST_Transform( ST_SetSRID(ST_MakePoint(8.8211, 47.2221), 4326), 3857)) "EPSG:3857 spherical mercator",
 +
    ST_AsText(ST_Transform( ST_SetSRID(ST_MakePoint(8.8211, 47.2221), 4326), 21781)) "EPSG:21781 CH1903/LV03"
 +
 +
=== Turn a latitude and longitude column into a geometry column ===
 +
 +
Often data is delivered where geometry is stored in separat columns, like in [[CSV]] format.
 +
 +
* First, make sure the latitude_column and longitude_column are of type NUMBER.
 +
* Then, make sure there is a column called the_geom which is of type POINT (and in the target [[CRS]], typically 4326)
 +
 +
Then call this SQL command which takes the latitude/longitude columns and creates a proper point geometry for each row in the database:
 +
  UPDATE table_name
 +
  SET the_geom = ST_SetSRID(ST_MakePoint(longitude_column, latitude_column),4326)
 +
 +
=== Join Lines ===
 +
 +
Line Join functionality. Takes connected separate lines which do no intersect and merges them into longer lines (Similar to FME LineJoiner, http://docs.safe.com/fme/2010/html/FME_Transformers/content/transformers/linejoiner.htm):
 +
 +
  SELECT (ST_Dump(ST_Linemerge(ST_Collect(geometry)))).geom from line;
 +
 +
Needs PostGIS >= 2.0. There's also the command "ST_CollectionHomogenize(<<geometry collection>>)" (needs PostGIS 2.x and GEOS >= 2.1.0).
 +
 +
=== Orientation of a Polygon ===
 +
[[Datei:Orientation_of_a_Polygon.jpg|200px|thumb|right|Orientation of a Polygon]]
 +
The main orientation of a polygon can be determined by the maximum distance of nodes of a polygon. To calculate this one can use ST_LongestLine() with the same geometry twice as input:
 +
 +
  SELECT ST_AsText( ST_Transform(ST_LongestLine(geom,geom),3857) ) as longestline
 +
  FROM ST_GeomFromText('POLYGON ((8.816893100738525 47.22353933792977, 8.817365169525146 47.223779787911226,
 +
  8.818191289901733 47.22307300758443, 8.818459510803223 47.223269741199125, 8.81878137588501 47.22301471600286,
 +
  8.817751407623291 47.222584813609146, 8.81751537322998 47.222796122001014, 8.817965984344482 47.22298557018802,
 +
  8.816893100738525 47.22353933792977))',4326) as geom;
 +
  -- Output: LINESTRING(981492.050352938 5978638.01719605,981702.252180723 5978552.02553924)
 +
 +
=== Clip all features from within a BBOX ===
 +
 +
Clip all features from within a BBOX (= Campus HSR Rapperswil; runs also in PostGIS Terminal):
 +
 +
  with bbox as (
 +
    select ST_Transform(ST_MakeEnvelope(8.813095,47.222229,8.819543,47.224415, 4326), 900913) AS way
 +
  ),
 +
  polygons as (
 +
    select ST_AsText(ST_Intersection(my.way,bbox.way)) as geom, osm_id, name, amenity
 +
    from osm_polygon my, bbox
 +
    where ST_Intersects(my.way,bbox.way)
 +
  ),
 +
  linestrings as (
 +
    select ST_AsText(ST_Intersection(my.way,bbox.way)) as geom, osm_id, name, amenity
 +
    from osm_line my, bbox
 +
    where ST_Intersects(my.way,bbox.way)
 +
  ),
 +
  points as (
 +
    select ST_AsText(my.way) as geom, osm_id, name,amenity
 +
    from osm_point my, bbox
 +
    where ST_Intersects(my.way,bbox.way)
 +
  )
 +
  select * from polygons
 +
  union
 +
  select * from linestrings
 +
  union
 +
  select * from points
 +
 +
=== Safe ST_UNION ===
 +
 +
=> https://trac.osgeo.org/postgis/attachment/ticket/4182/fun_geo_safe_union4.sql
  
 
[[Kategorie:PostGIS]] [[Kategorie:HowTo]]
 
[[Kategorie:PostGIS]] [[Kategorie:HowTo]]

Aktuelle Version vom 26. April 2022, 21:30 Uhr

See also:


Inhaltsverzeichnis

PostGIS-Datenbank mit Geometrie-Tabellen erzeugen

Quelle: BostonGIS.

Erläuterungen: "shell>" ist der Prompt einer DOS- oder Linux-Shell. "gisdb=#" ist der Prompt von psql, eingeloggt als gisdb-User.

PostgreSQL und PostGIS-Versionen kontrollieren
  • PostgreSQL-Version: SELECT version();
  • PostGIS-Version:
 SELECT postgis_full_version(); -- (bei aktiver Datenbank)
  • sonst template kontrollieren?
PostGIS-Extension installieren
 gisb=# CREATE EXTENSION postgis;
PostGIS-Tabelle erstellen (mit GEOMETRY-Typ und "type modifier (typmod)))
 CREATE TABLE geofoo (
   id int4 PRIMARY KEY, 
   name text,
   geom geometry('POINT', 21781)
 )
PostGIS-Index erzeugen und Installation überprüfen

PostGIS-Index erzeugen:

 shell> psql -U gisdb
 gisdb=# CREATE INDEX idx_towns_the_geom ON towns USING gist(the_geom);
 gisdb=# CREATE INDEX idx_towns_town ON towns USING btree(town);

Shapefiles in PostGIS importieren (shp2pgsql)

Man beachte beim Import von Shapefiles das .prj File. Ohne weitere Angaben wird der SRID auf -1 gesetzt. Siehe SRID (Quelle).

Shapefile-Daten in Datenbank geo1 laden:

 % shp2pgsql -s 21781 -I -W ISO-8859-1 gemeinden public.gemeinden > gemeinden.sql
 % psql -d geo1 -U geo1 -f gemeinden.sql

Eine Alternative zu shp2pgsql ist ogr2ogr von OGR.

PostGIS-Daten laden

  • Mittels (mit PostGIS ausgeliefertem) Kommandozeilen-Tool 'shp2pgsql'.
  • Mittels INSERT und WKT oder Konstruktoren (vgl. Snippet "Creating Geometry (...) Types" unten).

Creating Geometry and Geography Types

Allowed and non-functioning constructors for creating geometry type POINT (lon lat, +/-180 +/-90):

 /*0*/ SELECT ST_AsEWKT(ST_GeomFromText('POINT(-71.06 42.28)'))           -- Preferred simplest text form without SRID
 /*1*/ SELECT ST_AsEWKT(ST_GeomFromText('POINT(-71.06 42.28)', 4326))     -- Preferred for text form with SRID
 /*2*/ SELECT ST_AsEWKT(ST_GeomFromText('SRID=4326;POINT(-71.06 42.28)')) -- Note: srid 'inline'
 /*3*/ SELECT ST_AsEWKT(ST_GeomFromEWKT('SRID=4326;POINT(-71.06 42.28)')) -- Alternative to ST_GeomFromText()
 /*4*/ SELECT ST_AsEWKT('SRID=4326;POINT(-71.06 42.28)'::geometry)        -- With cast; prefer *1*/ST_GeomFromText
 /*5*/ SELECT ST_AsEWKT(ST_SetSRID('POINT(-71.06 42.28)'::geometry,4326)) -- With cast; prefer *1*/ST_GeomFromText
 /*6*/ SELECT ST_AsEWKT(ST_SetSRID(ST_MakePoint(-71.06, 42.28),4326))     -- Preferred symbolic form with EWKT

Allowed and non-functioning constructors for creating geography types (Varianten für geography types):

 /*1*/  -- SELECT ST_AsEWKT(ST_GeogFromText('POINT(-71.06 42.28)', 4326))       -- ERROR: unknown ST_GeogFromText()
 /*2a*/ -- SELECT ST_AsEWKT(ST_GeogFromText('SRID=4326;POINT(-71.06 42.28)'))   -- ERROR: no ST_AsEWKT: why not?
 /*2b*/ SELECT ST_AsText(ST_GeogFromText('SRID=4326;POINT(-71.06 42.28)'))      -- Preferred with srid 'inline'; (Hint: ST_AsEWKT doesn't work)
 /*3*/  -- SELECT ST_AsText(ST_GeogFromEWKT('SRID=4326;POINT(-71.06 42.28)'))   -- ERROR: no ST_GeogFromEWKT()
 /*4*/  SELECT ST_AsText('SRID=4326;POINT(-71.06 42.28)'::geography)            -- With cast; prefer *1*; (Hint: ST_AsEWKT doesn't work)
 /*5*/  --- SELECT ST_AsText(ST_SetSRID('POINT(-71.06 42.28)'::geography,4326)) -- ERROR: ok SetSRID w. geography unnecessary
 /*6*/  SELECT ST_MakePoint(-71.06, 42.28)                                      -- for geography types missing

Create Polygon given Bounding Box (BBox) Coordinates

Try this in http://labs.geometa.info/postgisterminal :

 SELECT ST_AsText(ST_Transform(ST_MakeEnvelope(8.795611, 46.872886, 9.674135, 47.675419, 4326), 3857)) AS geom, '#' AS label
 SELECT ST_AsText(ST_Transform(ST_SetSRID(ST_Envelope('LINESTRING(8.795611 46.872886, 9.674135 47.675419)'::geometry),4326), 3857)) AS geom, '#' AS label
 SELECT ST_AsText(ST_Transform(ST_SetSRID('BOX(8.795611 46.872886, 9.674135 47.675419)'::box2d, 4326), 3857)) AS geom, '#' AS label
 SELECT ST_AsText(ST_Transform(ST_SetSRID('BOX3D(8.795611 46.872886, 9.674135 47.675419)'::box3d, 4326), 3857)) AS geom, '#' AS label
 SELECT ST_AsText(ST_SetSRID(ST_Extent(way), 4326)) AS geom FROM osm_point

Curved Geometries

Curved geometries in PostGIS with notes for QGIS (and OGR).

Introduction: See GeoPackage Spec. (http://www.geopackage.org/spec110/#sfsql_intro) und PostGIS 2.5 Manual (https://postgis.net/docs/manual-2.5/. Citing PostGIS 2.5 Manual "4.1.3 SQL-MM Part 3": The SQL Multimedia Applications Spatial specification extends the simple features for SQL spec by defining a number of circularly interpolated curves. The SQL-MM definitions include 3DM, 3DZ and 4D coordinates, but do not allow the embedding of SRID information. The Well-Known Text extensions are not yet fully supported.". This diagram is interestingly showing the SQL-MM geometry types: http://www.dpi.inpe.br/terralib5/wiki/doku.php?id=wiki:documentation:devguide:geometry_module :

Iso sql mm geometry type hierarchy.png

These are the curved types (from the PostGIS 2.5 Manual):

  • CIRCULARSTRING is the basic curve type, similar to a LINESTRING. A single segment has second point as point on the arc (there are also closed circles, chained arcs)
  • COMPOUNDCURVE is a single, continuous curve that has both curved (circular) segments and linear segments.
  • CURVEPOLYGON is just like a polygon. The difference is that a outer/inner ring can take the form of a circular string, linear string or compound string.
  • MULTICURVE is a collection of curves, which can include linear strings, circular strings or compound strings.
  • MULTISURFACE is a collection of surfaces, which can be (linear) polygons or curve polygons.

Note that MultiCurve, MultiSurface are called abstract (and reserved) in the OGC model. Astonishingly in PostGIS, GeoPackage and QGIS it's used anyway.

QGIS: Supports these (CircularString, CompoundCurve, CurvePolygon, MultiCurve, MultiSurface). See https://www.qgis.ch/de/ressourcen/anwendertreffen/2015/curved-geometries-interlis-with-ogr-and-qgis .

OGR: Some functions to convert curves/arcs to linstrings.

Examples: (Note one could omit LINESTRING in 'LINESTRING(1 0, 0 1))' and also POLYGON here):

   SELECT ST_AsText('CIRCULARSTRING(0 0, 1 1, 1 0)'::geometry); 
   SELECT ST_AsText('COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 1 0),LINESTRING(1 0, 0 1))'::geometry);
   SELECT ST_AsText('CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0,2 0, 2 1,2 3,4 3),LINESTRING(4 3, 4 5, 1 4, 0 0)), CIRCULARSTRING(1.7 1, 1.4 0.4, 1.6 0.4, 1.6 0.5, 1.7 1))'::geometry);
   SELECT ST_AsText('MULTICURVE((0 0, 5 5),CIRCULARSTRING(4 0, 4 4, 8 4))'::geometry);
   SELECT ST_AsText('MULTISURFACE(CURVEPOLYGON(CIRCULARSTRING(0 0, 4 0, 4 4, 0 4, 0 0),LINESTRING(1 1, 3 3, 3 1, 1 1)),POLYGON((10 10, 14 12, 11 10, 10 10),(11 11, 11.5 11, 11 11.5, 11 11)))'::geometry);


Functions:

  • ST_IsCollection returns TRUE if the geometry type of the argument is either a GEOMETRYCOLLECTION, MULTI{POINT,POLYGON,LINESTRING,CURVE,SURFACE}, or a COMPOUNDCURVE. NOTE: This PostGIS-Doc. should also mention CURVEPOLYGON!
  • ST_HasArc: e.g. show if it's an art aber breaking up a compound curve into its segments
   SELECT ST_AsEWKT(a.geom), ST_HasArc(a.geom)
   FROM (SELECT (ST_Dump(p_geom)).geom AS geom
     FROM (SELECT ST_GeomFromEWKT('COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 1 0),(1 0, 0 1))') AS p_geom) AS b
   ) AS a;
   -- "CIRCULARSTRING(0 0,1 1,1 0)"    true
   -- "LINESTRING(1 0,0 1)"            false
  • ST_ForceCurve — Upcast a geometry into its curved type, if applicable.
   SELECT ST_AsText(
       ST_ForceCurve(
         'POLYGON((0 0, 5 0, 0 5, 0 0),(1 1, 1 3, 3 1, 1 1))'::geometry
       )
   );
   -- "CURVEPOLYGON((0 0,5 0,0 5,0 0),(1 1,1 3,3 1,1 1))"
  • ST_LineToCurve — Converts a LINESTRING/POLYGON to a CIRCULARSTRING, CURVEPOLYGON. If is not curved enough to clearly represent a curve the function will return the same input geometry.
  • ST_CurveToLine(geometry curveGeom, float tolerance, integer tolerance_type, integer flags) — Converts a CIRCULARSTRING/CURVEPOLYGON/MULTISURFACE to a LINESTRING/POLYGON/MULTIPOLYGON


PostGIS-Daten exportieren

Export a table with a geometry attribute to GeoJSON

Sources:

Query:

 :: File rowjsonb_to_geojson.sql ::
 /* ------------------------------------------------------------------
 Usage :
   SELECT rowjsonb_to_geojson(to_jsonb(mytable.*)) FROM mytable; 
 Parameters:
   JSONB - jsonb data type (aka one row)
   geom_column - name of geometry column (default 'geom')
 Credits: P. Ramsey, http://blog.cleverelephant.ca/2019/03/geojson.html
 ------------------------------------------------------------------ */
 CREATE OR REPLACE FUNCTION rowjsonb_to_geojson(
   rowjsonb JSONB, 
   geom_column TEXT DEFAULT 'geom')
 RETURNS TEXT AS $$
 DECLARE 
   json_props jsonb;
   json_geom jsonb;
   json_type jsonb;
 BEGIN
   IF NOT rowjsonb ? geom_column THEN
      RAISE EXCEPTION 'geometry column % is missing', geom_column;
   END IF;
   json_geom := ST_AsGeoJSON((rowjsonb ->> geom_column)::geometry)::jsonb;
   json_geom := jsonb_build_object('geometry', json_geom);
   json_props := jsonb_build_object('properties', rowjsonb - geom_column);
   json_type := jsonb_build_object('type', 'Feature');
   return (json_type || json_geom || json_props)::text;
 END; 
 $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT;

(Query - older solution):

 :: File togeojson.sql ::
 /* ---------------------------
 Convert following query output to GeoJSON:
   SELECT 
     st_transform(geom,4326) AS geom, -- geometry
     gid as id, "name", "level", pop_class -- "properties"
   FROM orte
   LIMIT 10 
 --------------------------- */
 SELECT jsonb_build_object(
   'type',     'FeatureCollection',
   'features', jsonb_agg(feature)
 ) AS geojson
 FROM (
   SELECT jsonb_build_object(
     'type',       'Feature',
     --'id',         id, -- 'id' is ignored e.g. by geojson.io
     'geometry',   st_asgeojson(geom)::jsonb,
     'properties', to_jsonb(row) - 'geom' -- removing elements with '-' 
   ) AS feature
   FROM ( -- input_table (beware columns named like keywords)
     SELECT 
       st_transform(geom,4326) AS geom, 
       gid as id, "name", "level", pop_class 
     FROM orte 
     LIMIT 10
   ) AS row
 ) AS features; 

Call it e.g. using psql like this:

 % psql -U <username> -d <dbname> -f togeojson.sql -t -o out.geojson

Ev. beautify output using geojson.io or http://geojsonlint.com/

Export a table with a geometry attribute to KML

PostGIS kennt die Funktion ST_AsKML(), die eine Geometrie KML-konform ausgibt. Allerdings werden nur diejenigen KML-Elemente der Geometrie selbst ausgegeben. Kopf- und Fusszeile des KML-Dokuments muss man selber hinzufügen. Dazu erzeugen wir nun eine Funktion askmldoc, die ein gültiges KML-Dokument ausgibt:

 CREATE OR REPLACE FUNCTION askmldoc(name text, description text, the_geom geometry) RETURNS text AS $$
 DECLARE
   result text;
 BEGIN
   result := '<?xml version="1.0" encoding="UTF-8"?>' || E'\n' ||
                   '<kml xmlns="http://www.opengis.net/kml/2.2">' || E'\n' ||
                   '<Document>' || E'\n' ||
                   '<name>' || name || '</name>' || E'\n' ||
                   '<description>' || description || '</description>' || E'\n\n' ||
                   '<Style id="defaultStyle">' || E'\n' ||
                   '  <LineStyle>' || E'\n' ||
                   '    <color>ff00ff00</color>' || E'\n' ||
                   '    <width>1</width>' || E'\n' ||
                   '  </LineStyle>' || E'\n' ||
                   '  <PolyStyle>' || E'\n' ||
                   '    <color>5f00ff00</color>' || E'\n' ||
                   '  </PolyStyle>' || E'\n' ||
                   '</Style>' || E'\n\n' ||
                   '<Placemark>' || E'\n' ||
                   '<styleUrl>#defaultStyle</styleUrl>' || E'\n';
         result := result || ST_AsKML(the_geom) || E'\n';
         result := result || 
                   '</Placemark>' || E'\n\n' ||
                   '</Document>' || E'\n' ||
                   '</kml>';
 RETURN result;
 END;
 $$ LANGUAGE plpgsql;

Exportieren Sie damit die Geodaten in der Tabelle, z.B. mit folgendem Befehl auf der System-Kommandozeile (cmd):

% psql -A -t -d gisdb -c "SELECT askmldoc('MyTable', 'A Comment...', the_geom) FROM mytable WHERE gid=1;" -o mykmlfile.kml


PostGIS-Daten darstellen

Übersicht:

PostGIS und Google Earth

PostGIS über http-Tool mit Google Earth verknüpfen:

Von PostGIS direkt nach Google Earth:

  • Man starte psql (Beispiel mit Box um Victoria, BC, Kanada):
 -- Set output to unaligned
 \a
 -- Show only tuples
 \t
 -- Dump query to a file
 \o my_kml_file.kml
 -- Run your query
 SELECT askml('my_line', 'description', 'SRID=3005;LINESTRING( 1190000
 390000, 1200000 390000, 1200000 380000, 1190000 380000, 1190000 390000
 )'::geometry);
 -- Flush output your file and close the output stream
 \o
  • Dann Doppelklick auf die KML-Datei und Google Earth startet!
  • Hinweis: Das KML muss ev. mit KML-Kopf und -Fuss ergänzt werden.
  • Beispiel für einen Ausschnitt der Schweiz:
 SELECT AsKML('the_geom', 
              'SRID=21781;LINESTRING(480000 300000, 835000 300000, 835000 70000, 480000 70000, 480000 300000)'::geometry
        );
  • Beispiel für gisdb:
 SELECT ST_AsKML(geomunion(transform(the_geom,4326))) as the_geom from towns where town='BOSTON';

Spatial Relationships and ST_Relate

ST_Relate is based on Dimensionally Extended Nine-Intersection Model (DE-9IM) (or Clementini-Matrix). It covers all known spatial relationships. The most typical spatial relationsships (and it's opposites) got own functions, like:

  • ST_Within != ST_Contains (oder ST_Covers (*))
  • ST_Covers != ST_CoveredBy
  • ST_Intersects != ST_Disjoint

(*) Note: Prefer ST_Covers over ST_Contains if lines on boundaries count as „inside“ (Source: Martin Davis‘ blog post: http://lin-ear-th-inking.blogspot.ch/2007/06/subtleties-of-ogc-covers-spatial.html )


Documentation

ST_Relate(geom_a, geom_b, '<<pattern>>')

Examples of pattern values:

  • ST_Equals:
    • 0******** (for Point/Point)
    • T*F**FFF* (for any geom?; from Esri Webhelp)
  • ST_Contains:
    • T*****FF* (for any geom?; from Esri Webhelp)
  • ST_Crosses:
    • T*T****** (for Point/Line, Point/Area, and Line/Area situations)
    • T*****T** (for Line/Point, Area/Point, and Area/Line situations)
    • 0******** (for Line/Line situations)
    • 0F1FF0102 (for strict Line/Line situation like an 'X')
  • ST_Touches (From [5]:"The allowable DE-9IM Intersection Matrices for the two geometries are:" ???):
    • FT*******  ?
    • F**T*****  ?
    • F***T****  ?
  • ST_Contains(polygon_a, polygon_b) => ST_Relate(polygon_a.geom, table_b.geom, '2********')
  • ST_Intersects => 'T********' (or within polygon/polygon: 'T*F**F***' ?)
  • ST_Overlaps
    • '1********' (for Area/Area)
  • ST_Within
    • T*F**F*** (for any geom?; from Esri Webhelp)
  • disjoint line/point => 'FF*FF****'
  • Census blocks and voting district interiors that do not intersect "F" (interior/interior) but have a common linear "1" boundary (boundary/boundary) => (polygonA, polygonB) => 'F***1****'
  • What about this topology rule: geometries that neither intersect nor overlap. The only allowed shared points between geometries are boundary-boundary intersections. boundary-interior-boundary (touch interior) intersections are not allowed. (From: Rule.java, Jaspa project [6]

ST_Relate(geom_a, geom_b)

Haben Sie gewusst, dass es auch eine Variante ST_Relate(geomA, geomB) gibt, die die DE-9IM „berechnet“ gegeben einen Geometrie-Input? Hier z.B. zwei sich schneidende Linien (als X) als Input:

 SELECT ST_Relate(PA,PB)
 FROM (
   SELECT ST_GeomFromText('LINESTRING(1 3, 9 6)') AS PA, 
          ST_GeomFromText('LINESTRING(3 6, 7 3)') AS PB
   ) AS foo
 => 0F1FF0102

PostGIS und Rasterdaten

Ab PostGIS Version 2.0 (geplant 1. Quartal 2011) gibt es neu einen Raster-Datentyp.

Siehe PostGIS Raster.

PostGIS und Koordinatenreferenzsystem-Angaben

Der SRID (projection identifier) wird an drei Orten verwaltet:

  • In der geometry column: select SRID(the_geom) from mytable limit 5;
  • As a constraint on the table for the geometry column: \d mytable
  • In the geometry_columns system table: select * from geometry_columns;

Mit der Methode UpdateGeometrySRID([<schema_name>], <table_name>, <column_name>, <srid>) kann man die Kolonne aktualisieren.

Das wohl bekannteste Koordinatenreferensystem ist wohl WGS 84 (long/lat), das von GPS und KML geprägt ist und den Identifier 'EPSG:4326' hat.

Liste von SRID/SRS/CRS: http://www.spatialreference.org/

Beim Aufruf von GDAL/OGR wird die GDAL_DATA environment Variable benötigt.

Erkennen und Eliminieren von Sliver Polygonen und Spikes

Implementiert als sog. 'Stored Procedures'.

Siehe auch PostGIS-Beispiele.

Eliminate sliver polygons

Sources:


Given a polygon table that has many small areas and holes. How to remove "small" areas and holes (smaller than a given area in m2)?

Remarks:

 CREATE OR REPLACE FUNCTION Filter_Rings(geometry,float)
 RETURNS geometry AS
 $$
 SELECT ST_MakePolygon(c.outer_ring, d.inner_rings) as final_geom
  FROM (/* Get outer ring of polygon */
        SELECT ST_ExteriorRing(b.the_geom) as outer_ring
          FROM (SELECT (ST_DumpRings($1)).geom As the_geom, path(ST_DumpRings($1)) as path) b
          WHERE b.path[1] = 0 /* ie the outer ring */
        ) c,
       (/* Get all inner rings > a particular area */
        SELECT ST_Accum(ST_ExteriorRing(b.the_geom)) as inner_rings
          FROM (SELECT (ST_DumpRings($1)).geom As the_geom, path(ST_DumpRings($1)) as path) b
          WHERE b.path[1] > 0 /* ie not the outer ring */
            AND ST_Area(b.the_geom) > $2
        ) d
 $$
 LANGUAGE 'sql' IMMUTABLE;
 CREATE OR REPLACE FUNCTION Filter_Rings(geometry GEOMETRY, param FLOAT)
 RETURNS geometry AS
 $$
 SELECT ST_MakePolygon(c.outer_ring, d.inner_rings) as final_geom
  SELECT * FROM foo WHERE attr > param; -- statt $2 steht hier ein benanntes Fn.-Argument
  ...
 $$
 LANGUAGE 'sql';

Usage example:

 % SELECT ST_AsText( 
   Filter_Rings( 
     ST_PolyFromText(
       'POLYGON((10 10,10 20,20 20,20 10,10 10),(0 0,0 1,1 1,1 0,0 0),(5 5,5 7,7 7,7 5,5 5))'
     ) ,1::float
   )
 );
 % "POLYGON((10 10,10 20,20 20,20 10,10 10),(5 5,5 7,7 7,7 5,5 5))"

Shorter alternative to Filter_Rings:

 CREATE OR REPLACE FUNCTION Filter_Rings2(geometry,float) RETURNS geometry AS
 $$ SELECT ST_BuildArea(ST_Collect(a.geom)) as final_geom
         FROM ST_DumpRings($1) AS a
           WHERE a.path[1] = 0 OR
                 (a.path[1] > 0 AND ST_Area(a.geom) > $2)
 $$
 LANGUAGE 'sql' IMMUTABLE;

With following restrictions: Possibly slower and squashes 3D geometries to 2D. (Source: http://postgis.refractions.net/pipermail/postgis-users/2009-January/022325.html)

Clean topology

Remarks: Similar to ArcGIS' CLEAN. See also Eliminate sliver polygons.

See:

PostGIS optimieren

Siehe auch:

Version:

  • 1.4 and GEOS 3.1 will bring prepared geometries which will make things faster.

Fast map intersection by spatial partioning https://github.com/gbb/fast_map_intersection

Reorder on disk:

  • import the data (into a temporary table)
  • and then sort the data using some sort of spatial key. A centroid should be OK but consider using a larger bucket such as a gridded area.
 > CREATE TABLE gis_roads 
     AS SELECT * FROM tmp_gis_roads 
     ORDER BY CENTROID(the_geom);

Cluster on disk:

 > CREATE INDEX idx_mytable_the_ON mytable 
     USING GIST(the_geom);
 > CLUSTER mytable 
     USING idx_mytable_the_geom;

Don'ts:

  • Don't use just one big diagonal lines. Break it down because one diagonal line has a huge bounding box (which is what the index works on) while 20 equivalent shorter lines have a much smaller combined box area.
  • Don't use this construction: ST_Intersects(Buffer(GeomA, D), GeomB) Use this one: ST_DWithin(GeomA, GeomB, D)

Weitere Beispiele

Siehe auch:

Fragen und Antworten

Was ist der Unterschied von geometry und geography type?
Geography type besteht aus geografischen Koordinaten (z.B. lat/lon) und wird für grossflächige Daten verwendet, die z.B. "länderübergreifend" sind.
Wie wird ein Geometrie-Attribut erzeugt?
als separates Statement, gleich nach CREATE TABLE... Mit AddGeometryColumn. Ab PostGIS 1.5 ist es endlich möglich, den konkreten Geometrie-Typ direkt im CREATE STATEMENT anzugeben.
Was bedeutet SRID und wie kann man das ändern?
SRID bedeutet Spatial Reference ID und ist ein Fremdschlüssel zum Koordinatenreferenzsystem. SRID information is stored *both* in the geometry_columns table (srid column), and inside each geometry itself. Ist opbligatorischer Parameter u.a. von AddGeometryColumn oder beim Geometrie-Erzeugen. Kann mit UpdateGeometrySRID geändert werden.
Was bedeutet der Error "addgeometrycolumn(...) does not exist"?
Die PostGIS-Funktionen sind nicht geladen. Man installiere die PostGIS-Erweiterung (falls nicht schon geschehen) und erzeuge eine neue PostGIS-Datenbank (siehe "Wie erzeugt man eine neue PostGIS-Datenbank?").
Wie erzeugt man eine neue PostGIS-Datenbank?
1. Falls noch keine 'normale' (=Nicht-GIS-) Datenbank vorhanden ist, "create a new database and choose the template_postgis as your template". 2. Falls eine 'normale' Datenbank mit Tabellen etc. bereits existiert, "spatially enable your database by running the following two SQL files: a) Postgresql\8.2\share\contrib\lwpostgis.sql and b) Postgresql\8.2\share\contrib\spatial_ref_sys.sql. lwpostgis.sql benötigt plpgsql language support installed (Syntax: createlang plpgsql yourexistingdb).
Was ist das für eine Zahl, die ST_distance() zurückgibt?
Siehe PostGIS Wiki.
Was ist der Unterschied zwischen ST_Overlaps, ST_Crosses, ST_Intersects und ST_Contains?
Siehe PostGIS-Doku. zu ST_Overlaps.
Was ist der Unterschied zwischen ST_Extent und ST_Expand?
Siehe [8].
Was für Mass-Einheiten (units of measurement) benützen die PostGIS-Funktionen?
Siehe PostGIS Wiki.
Wie erhält man Distanz-Masse in Meter?
Siehe PostGIS Wiki.
Was sind die Grenzen der 3D-Typen und -Funktionen in PostGIS?
(1) 3D volumetric objects are not supported. (2) 3D non-volume are supported partly -- e.g. a 2D polygon in 3 d space, a line, point in 3d space. (3) Spatial relationships however only consider the spatial component of the x, y plane.
MySQL hat auch eine Spatial Extension und ist erst noch bei Internet-Service-Providern verbreiteter; was ist der Unterschied zu PostGIS?
Drei Gründe: (1) Operationen werden nur mittels Bounding-Box gemacht (hips-between-geometries.html, (2) es gibt keine Unterstützung für Koordinatensysteme, (3) die Spatial Extension ist praktisch eine One-Man-Show.

Snippets

Find centre and radius of a curve defined by three points

“COGO: Finding centre and radius of a curve defined by three points (PostGIS)” by Simon Greener: https://spatialdbadvisor.com/postgis_tips_tricks/278/postgis-cogo-finding-centre-and-radius-of-a-curve-defined-by-three-points

Get area of a polygon which has SRID EPSG:4326 (WGS84/GPS)

All area functions are squared by definition. But a sq degree varies in its area depending on latitude. Transform your data using a SRID referring to an equal-area (or near enough to equal area for your purposes) projection in the query for area:

 -- Get area from a geometry given the geometry is in 4326.
 % SELECT ST_Area2d(ST_transform(geometry, <SRID>));

Extent of Database

If you want to find out the EXTENT parameters of your database, run this command:

 # su postgres
 # psql -d <yourdatabase> -c "SELECT extent(geometry) FROM <yourtable>"

Replace <your...> with the relevant database and table names.

Nearest Neighbors

Find n (e.g. 5) nearest neighbors for given point coordinate (a classic query). From '[postgis-users] Find n Nearest Neighbors for given Point using PostGIS?'.

 CREATE TABLE geoname ( -- (from geonames.org) 
   geoname varchar(255),
   geonameid integer,
   geom geometry
 );
 -- PRIMARY KEY UNIQUE BTREE index for geonameid
 
 UPDATE geoname SET geom = ST_SetSRID(ST_Point(longitude,latitude) 4326);
 ...
 CREATE INDEX geom_index ON geoname USING GIST (geom);) 
 CLUSTER geom_index ON geoname;

Possible solution (inspired by [9]):

 SELECT 
   start.asciiname, 
   ende.asciiname, 
   distance_sphere(start.geom,ende.geom) as distance
 FROM geoname AS start, geoname AS ende 
 WHERE start.geonameid = 2950159 
   AND start.geonameid <> ende.geonameid 
   AND ST_DWithin(start.geom, ende.geom, 30) 
 ORDER BY distance LIMIT 5;

There will be a KNNGIS index probably in PostgreSQL 9.1 which makes such queries even faster. It returns the found tupels in decreasing order until LIMIT without computing and sorting whole set before returning it.

The way currently is to write a stored procedure that expands the search if it fails to get the number of results. So make a set returning function with the body something like this (pseudo code) and you should get good performance:

 radius := 0.01; -- assuming degrees
 loop
  select into cnt count(*) from mytable
   where expand(mypnt, radius) && the_geom limit 5;
  if found and cnt = 5 or radius > maxradius then
    for rr in select * from mytable
                where expand(mypnt, radius) && the_geom limit 5
      loop
        return rr;
      end loop;
    return;
  else
    radius := radius * 2;
  end if;
 end loop;

Find nearest existing point on a linestring given a point (Snapping 1)

A user clicks somewhere on a polyline and we want to get the closest existing point that is part of the representation of the line. You can replace 0.01 with some other tolerance number to compensate for the fact that a person probably won't click right on the line.

Ab PostGIS 1.3.5 (wegen generate_series() )?

 SELECT P.gid, P.geom As pclosest
 FROM (SELECT 
   gid, 
   ST_PointN(the_geom, generate_series(1, ST_NPoints(the_geom))) AS geom
 FROM point_table) AS P
 WHERE ST_Dwithin(P.geom, ST_SetSRID(ST_MakePoint(lon,lat),somesrid), 0.01)

Find nearest interpolated point on a linestring given a point (Snapping 2)

See http://blog.cleverelephant.ca/2008/04/snapping-points-in-postgis.html

Find nearest (whole) linestring given a point (Snapping 3)

Basic principle given mypos => 'SRID=3005;POINT(1011102 450541)'

Solution using distance within max. radius e.g. in meters (ST_DWithin):

 SELECT geom, name, etc FROM mytable
 WHERE ST_DWithin(mypos, geom, 50) -- max. distance 50m
 ORDER BY ST_Distance(mypos, geom)
 LIMIT 1

and/or solution using 50 nearest objects (kNN index). Note that's returng 50 objects even if they are far away.

 WITH tmp AS (
   SELECT geom, name, ST_Distance(geom, mypos) AS distance
   FROM mytable
   ORDER BY geom <#> mypos -- prefering <#> over <-> when lines are involved.
   LIMIT 50
 )
 SELECT * FROM tmp 
 ORDER BY distance 
 LIMIT 1

See also the examples at PostGIS Terminal Examples.


Find and add nearest hydrant to a parcel

Given parcels (or other kind of (point?) feature class) and hydrants (or other kind of (point?) feature class) as input (e.g. imported from external files), find and add closest hydrant to parcel. Output: All parcel attributes incl. geometry, hydrant_id and distance from hydrant to parcel.

This is an example of a lateral join driving a nearest-neighbor distance calculation. Needs PostgreSQL 9.5 / PostGIS 2.x. (Credits: https://twitter.com/pwramsey/status/686980049967443968 )

 SELECT
   -- Take all attributes from parcels incl. pid and geom:
   parcels.*,
   -- Keep the hydrant id around, might be useful later:
   hydrants.id AS hydrant_id,
   -- Calculate distance over the spheroid using geography distance:
   ST_Distance(geography(hydrants.geom), geography(parcels.geom)) AS distance
 FROM 
   -- For this data, removing the duplicate parcel geometries and
   -- null attributed parcels cuts the size by 2/3 so a good
   -- performance improvement to filter records here:
   (SELECT DISTINCT ON (geom) *
    FROM parcels_imported 
    WHERE pid IS NOT NULL) AS parcels 
   -- Cross join takes every combination of join records
   -- without restriction, which is fine since we're limiting
   -- the hydrants side of the join to just one return record
   -- per parcel:
   CROSS JOIN LATERAL
     -- The guts of the query, actually do a KNN query for each
     -- candidate parcel on the other side of the join:
     (SELECT id, geom
      FROM hydrants_imported
      ORDER BY parcels.geom <-> geom -- take here cartesian CRS like webmercator 
      LIMIT 1) AS hydrants 
   ON TRUE;

PostGIS views

Q: I like to polygonize linestrings in a view to see the result in qgis. But qgis need gid. How can I create the gid column after an aggregate function in a view?

A1. In postgresql 8.3 you can force oid creation if WITH OIDS is specified when the table is created:

 CREATE VIEW myview AS SELECT oid AS gid, the_geom FROM mytable;

A2. Create a sequence and use then nextval(sequence_name) to get a primary key.

A3. Grab the key from the underlying geometry table, & use it for a key column in your aggregate view, eg. "create view v1 as select min(key) as key" ...

A4. Use window function "SELECT row_number() OVER() As generated_id, ... FROM ..."

Standard direction for polylines

Q: Is there a standard direction for polylines?

A: There isn't a standard, but you can force your polygons into an orientation with ST_ForceRHR, which forces a clockwise orientation for exterior rings and anti-clockwise for interior rings.

Portion of river geometry along administrative boundary

Q: I have a rivers linestring table and a political boundaries (counties) table polygon. I am wishing to extract the portion of the river geometry that lies along a county border. Of course, these two layers don't have the same individual point precision, so I suspect some buffering is necessary.

A: See following SQL query:

 intersection(
      buffer(exteriorring(geometryn(n.geom,1)),0.02),
      r.geom
 )

So turn the county polygon into an exterior ring and then buffer it out a bit. Then see what portions of the river intersect that buffered polygon.

Convert multipolygon to polygon

Q: How to convert multipolygon to polygon geometries

A: If multipolygon only have one member, then st_geometryn(geom, 1) will do it. Otherwise, look at function st_dump().

Getting all polygon's coordinates

Q. Given a geometry attribute with polygons, how can one get all the polygon's coordinates (as a tupel set)? E.g. either

 1 POINT(1 2)
 2 POINT(2 2)
 3 POINT(3 2)
 ...

... or even better:

     X    Y
 1  1     2
 2  2     2
 3  3     2
 ...
 

1. Maybe a "SELECT ST_AsText(polygon_column) FROM Table" is enough

Output: POLYGON( (1 1, 2 2, 3 3) )

2.

 SELECT 
   DISTINCT ST_X(ST_POINTN(foo.lines, foo.index)) AS x,
   ST_Y(ST_POINTN(foo.lines, foo.index)) AS y 
 FROM (SELECT thelines.lines, generate_series(1,thelines.n) AS index 
         FROM (SELECT polys.geo AS lines, ST_NUMPOINTS(polys.geo) AS n 
                 FROM (SELECT ST_BOUNDARY((ST_DUMP( the_geom )).geom) AS geo 
                         FROM "YOURTABLE" WHERE YOURQUERY
                     ) AS polys
             ) AS thelines
      ) AS foo;
 SELECT DISTINCT 
 -- foo.index,
 ST_X(ST_POINTN(foo.lines, foo.index)) AS x,
 ST_Y(ST_POINTN(foo.lines, foo.index)) AS y 
 FROM (SELECT thelines.lines, generate_series(1,thelines.n) AS index 
         FROM (SELECT polys.geo AS lines, ST_NUMPOINTS(polys.geo) AS n 
                 FROM (SELECT ST_BOUNDARY((ST_DUMP( the_geom )).geom) AS geo 
                         FROM "gemeinden" WHERE gid=4 
                     ) AS polys
             ) AS thelines
      ) AS foo

Should work for multipolygons and polygons.


3. You can import a plpgsql implementation of this until it's implemented in C (http://trac.osgeo.org/postgis/attachment/ticket/76/my_st_dump_points.2.sql). [[Media:]]

Usage:

 SELECT ST_AsText( (My_ST_DumpPoints( poly )).geom )
 FROM (
   SELECT 'POLYGON((0 0, 1 1, 1 0, 0 0))'::geometry AS poly
 ) AS foo;
 ST_AsText ------------
 POINT(0 0)
 POINT(1 1)
 POINT(1 0)
 POINT(0 0)
 (4 rows)

Split a polygon to N linestrings

Source: http://postgis.refractions.net/pipermail/postgis-users/2010-January/025818.html

 CREATE TEMP TABLE mypolygontable AS
 SELECT
 'MULTIPOLYGON(((0 0, 0 4, 4 4, 4 0, 0 0),
                (1 1, 1 2, 2 2, 2 1, 1 1)),
               ((5 5, 5 6, 6 6, 6 5, 5 5))
              )'::geometry geom;
 -- make line segments from every startpoint and endpoint
 SELECT ST_AsText( ST_MakeLine(sp,ep) )
 FROM
   -- extract the endpoints for every 2-point line segment for each
 (SELECT
    ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
    ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
  FROM
     -- extract the individual linestrings
    (SELECT (ST_Dump(ST_Boundary(geom))).geom
     FROM polygons
     ) AS linestrings
 ) AS segments;
     st_astext
 ---------------------
 LINESTRING(0 0,0 4)
 LINESTRING(0 4,4 4)
 LINESTRING(4 4,4 0)
 LINESTRING(4 0,0 0)
 LINESTRING(1 1,1 2)
 LINESTRING(1 2,2 2)
 LINESTRING(2 2,2 1)
 LINESTRING(2 1,1 1)
 LINESTRING(5 5,5 6)
 LINESTRING(5 6,6 6)
 LINESTRING(6 6,6 5)
 LINESTRING(6 5,5 5)
(12 rows)


More general solution:

Aggregate linestrings

Q. Trying to aggregate linestrings together based on their attributes and the fact that they are touching each other.

A. There is a ST_Linemerge function that takes a collection of linestrings as an argument and merges them together into a multilinestring.

With a table "thetable" having an "attrib" field and a "geom" field, one could write the query like this :

 SELECT attrib, (st_dump(merged_geom)).geom
 FROM (
   SELECT attrib, ST_Linemerge(ST_Collect(geom)) AS merged_geom
   FROM thetable
   GROUP BY attrib
 ) AS subq;

Polygons that don't overlap

Q. I have a table with a MULTIPOLYGON field. Some of these records have self intersections and other problems that make ST_IsValid return false. The polygons overlap a lot and I wanted to generate another table that will be the union of all polygons. The table looks like: id, the_geom. What would be the best way to union all of the geometries into a new table where there is only POLYGONS that do not overlap?

A. See SQL:

 SELECT id, ST_Union(ST_Accum(CASE WHEN ST_IsValid(a.the_geom) = 't' THEN a.the_geom ELSE st_buffer(a.the_geom,0) end) AS the_geom
 FROM (SELECT id, the_geom
   FROM mytable
   WHERE the_geom IS NOT NULL
     AND geometrytype(the_geom) != 'GEOMETRYCOLLECTION'
   UNION ALL
      SELECT gid, (ST_Dump(the_geom)).geom AS the_geom
        FROM mytable
        WHERE the_geom IS NOT NULL
          AND geometrytype(the_geom) = 'GEOMETRYCOLLECTION'
      ) AS a;

The split based on GEOMETRYCOLLECTION is that when I first tested the STR indexed ST_Union I found that performance was affected if the union set includes an mpoly.

Transform coordinates into other reference systems

The following statement "abuses" PostGIS to simply transforms the given point coordinates into other reference systems and prints the result to the console. This is useful for example to double-check self coded transformation algorithms.

 SELECT
   ST_AsText(ST_Transform( ST_SetSRID(ST_MakePoint(8.8211, 47.2221), 4326), 4326)) "EPSG:4326",
   ST_AsText(ST_Transform( ST_SetSRID(ST_MakePoint(8.8211, 47.2221), 4326), 3857)) "EPSG:3857 spherical mercator",
   ST_AsText(ST_Transform( ST_SetSRID(ST_MakePoint(8.8211, 47.2221), 4326), 21781)) "EPSG:21781 CH1903/LV03"

Turn a latitude and longitude column into a geometry column

Often data is delivered where geometry is stored in separat columns, like in CSV format.

  • First, make sure the latitude_column and longitude_column are of type NUMBER.
  • Then, make sure there is a column called the_geom which is of type POINT (and in the target CRS, typically 4326)

Then call this SQL command which takes the latitude/longitude columns and creates a proper point geometry for each row in the database:

 UPDATE table_name 
 SET the_geom = ST_SetSRID(ST_MakePoint(longitude_column, latitude_column),4326)

Join Lines

Line Join functionality. Takes connected separate lines which do no intersect and merges them into longer lines (Similar to FME LineJoiner, http://docs.safe.com/fme/2010/html/FME_Transformers/content/transformers/linejoiner.htm):

 SELECT (ST_Dump(ST_Linemerge(ST_Collect(geometry)))).geom from line;

Needs PostGIS >= 2.0. There's also the command "ST_CollectionHomogenize(<<geometry collection>>)" (needs PostGIS 2.x and GEOS >= 2.1.0).

Orientation of a Polygon

Orientation of a Polygon

The main orientation of a polygon can be determined by the maximum distance of nodes of a polygon. To calculate this one can use ST_LongestLine() with the same geometry twice as input:

 SELECT ST_AsText( ST_Transform(ST_LongestLine(geom,geom),3857) ) as longestline
 FROM ST_GeomFromText('POLYGON ((8.816893100738525 47.22353933792977, 8.817365169525146 47.223779787911226, 
 8.818191289901733 47.22307300758443, 8.818459510803223 47.223269741199125, 8.81878137588501 47.22301471600286, 
 8.817751407623291 47.222584813609146, 8.81751537322998 47.222796122001014, 8.817965984344482 47.22298557018802, 
 8.816893100738525 47.22353933792977))',4326) as geom;
 -- Output: LINESTRING(981492.050352938 5978638.01719605,981702.252180723 5978552.02553924)

Clip all features from within a BBOX

Clip all features from within a BBOX (= Campus HSR Rapperswil; runs also in PostGIS Terminal):

 with bbox as (
   select ST_Transform(ST_MakeEnvelope(8.813095,47.222229,8.819543,47.224415, 4326), 900913) AS way
 ),
 polygons as (
   select ST_AsText(ST_Intersection(my.way,bbox.way)) as geom, osm_id, name, amenity
   from osm_polygon my, bbox 
   where ST_Intersects(my.way,bbox.way)
 ),
 linestrings as (
   select ST_AsText(ST_Intersection(my.way,bbox.way)) as geom, osm_id, name, amenity
   from osm_line my, bbox 
   where ST_Intersects(my.way,bbox.way)
 ),
 points as (
   select ST_AsText(my.way) as geom, osm_id, name,amenity
   from osm_point my, bbox 
   where ST_Intersects(my.way,bbox.way)
 )
 select * from polygons
 union
 select * from linestrings
 union
 select * from points

Safe ST_UNION

=> https://trac.osgeo.org/postgis/attachment/ticket/4182/fun_geo_safe_union4.sql