PostGIS Terminal Examples: Unterschied zwischen den Versionen

Aus Geoinformation HSR
Wechseln zu: Navigation, Suche
(Filter multivalued tags from OpenStreetMap)
K (Adresse eines Shops)
(5 dazwischenliegende Versionen desselben Benutzers werden nicht angezeigt)
Zeile 179: Zeile 179:
 
=== Alle 4000er Berggipfel der Schweiz ===
 
=== Alle 4000er Berggipfel der Schweiz ===
 
  <nowiki>
 
  <nowiki>
 +
  CREATE OR REPLACE FUNCTION as_numeric(text) RETURNS NUMERIC AS $$
 +
  -- Inspired by http://stackoverflow.com/questions/16195986/isnumeric-with-postgresql/16206123#16206123
 +
  DECLARE test NUMERIC;
 +
  BEGIN
 +
      test = $1::NUMERIC;
 +
      RETURN test;
 +
  EXCEPTION WHEN others THEN
 +
      RETURN -1;
 +
  END;
 +
  $$ STRICT
 +
  LANGUAGE plpgsql IMMUTABLE;
 +
 +
  SELECT ST_AsText(way) AS geom, osm_id||', '||coalesce(name,'No name')||', '||as_numeric(ele)::int||'m' AS label
 +
  FROM osm_point
 +
  WHERE tags @> 'natural=>peak' AND as_numeric(ele) between 4000 and 9999
 +
  ORDER BY as_numeric(ele) desc;
 +
 +
  oder:
 +
 
   SELECT ST_AsText(way) AS geom, name||','||ele AS label
 
   SELECT ST_AsText(way) AS geom, name||','||ele AS label
 
   FROM osm_point
 
   FROM osm_point
Zeile 377: Zeile 396:
 
   AND tags ? 'wikipedia'
 
   AND tags ? 'wikipedia'
 
   ORDER BY way <-> start.geom
 
   ORDER BY way <-> start.geom
   LIMIT 15 </nowiki>
+
   LIMIT 15; </nowiki>
 +
 
 +
 
 +
=== Adresse eines Shops === 
 +
* Adresse eines Shops mit Attributen von seinen Umliegenden Objekten, wenn er selber keines hat.
 +
* Hier am Beispiel von Bankomaten (ATM):
 +
  <nowiki> select
 +
    atm.osm_id,
 +
    atm.name,
 +
    atm.tags,
 +
    coalesce((atm.tags->'addr:street'), building.street, anyobject.street) as street,
 +
    coalesce((atm.tags->'addr:housenumber'), building.housenumber, anyobject.housenumber) as housenumber,
 +
    st_transform(atm.way,4326) as geom
 +
  from osm_point as atm
 +
  left join lateral (
 +
    select (tags->'addr:street') as street, (tags->'addr:housenumber') as housenumber
 +
    from osm_polygon
 +
    where (tags->'building') is not null and (tags->'addr:street') is not null
 +
    and st_within(atm.way, way)
 +
  ) as building on true
 +
  left join lateral (
 +
    select (tags->'addr:street') as street, (tags->'addr:housenumber') as housenumber
 +
    from osm_polygon
 +
    where (tags->'addr:street') is not null
 +
    order by atm.way <-> way
 +
    limit 1
 +
  ) as anyobject on true
 +
  where tags @> 'amenity=>atm'::hstore
 +
  order by street; </nowiki>
  
 
== Special SQL Queries ==
 
== Special SQL Queries ==
Zeile 499: Zeile 546:
 
   */</nowiki>
 
   */</nowiki>
  
This is some text code showing step-by-step how the multivalue string is being converted to an array:
+
This is some test code showing step-by-step how the multivalue string is being converted to an array:
 
   <nowiki>with osm_poi_raw (id, name, cuisine_str) as ( values  
 
   <nowiki>with osm_poi_raw (id, name, cuisine_str) as ( values  
 
     (1, 'Reschti 1', 'italian;swiss;pizzeria'),  
 
     (1, 'Reschti 1', 'italian;swiss;pizzeria'),  
Zeile 582: Zeile 629:
  
 
== Statistics ==  
 
== Statistics ==  
 +
 +
=== Youngest node in the database ===
 +
 +
  select
 +
    osm_id,
 +
    name,
 +
    (tags->'osm_timestamp')::timestamptz as timestamp,
 +
    (tags->'osm_version')::int as version,
 +
    (tags->'osm_uid')::int as uid,
 +
    (tags->'osm_changeset')::int as changeset,
 +
    way
 +
  from osm_point
 +
  where osm_id=(select max(osm_id) from osm_point)
 +
  
 
=== Keys containing problematic colons ===
 
=== Keys containing problematic colons ===

Version vom 23. September 2020, 01:47 Uhr

Here are some examples queries for the PostGIS Terminal. All SQL queries work in the terminal (like this) but mostly also on any PostGIS database created with osm2pgsql.

See also PostGIS - Tipps und Tricks.

Fun with SQL Queries

Checking if a road is mapped twice

Figuring out if a road is mapped twice with the same nodes. So find all line objects that share at least one node. Use following query while zoomed in at something like zoom level 18 or 19:

with tmp as (
  select osm_id, way
  from osm_line 
  where st_contains(mapextent(), way)
  and osm_id > 0
)
select 
  st_astext(a.way) as geom, 
  a.osm_id||' '||b.osm_id as label
from 
  tmp as a,
  tmp as b
where st_overlaps(a.way, b.way);

Find nearest linestrings with key 'highway' given a point

WITH mypos AS (
  SELECT ST_SetSRID(ST_MakePoint(_mouse_x, _mouse_y), 3857) AS geom
),
tmp AS (
  SELECT osm_id, way, name, ST_Distance(way, (SELECT geom FROM mypos)) AS distance
  FROM osm_line
  WHERE (tags->'highway') IS NOT NULL
  ORDER BY way <#> (SELECT geom FROM mypos)
  LIMIT 100
)
SELECT 
  ST_AsText(way) geom,
  COALESCE(name, '')||' '||osm_id AS label
FROM tmp 
ORDER BY distance 
LIMIT 4

All Roundabouts

-- This query works but all attributes from joined lines in first query are "lost".
-- NOTE: Add "ST_Within(way, (SELECT way FROM osm_polygon WHERE osm_id=-51701))" 
-- for exact results within Switzerland but put the subquery out of ST_Within into the 
-- from-clause instead since the presence of a subquery is blocking the inlining of ST_Within.
WITH lines AS (
    SELECT 
        osm_id as id, 
        way as geom 
    FROM osm_line
    WHERE 
      (tags @> hstore('junction','roundabout') 
      OR tags @> hstore('highway','mini_roundabout'))
      AND NOT (tags ? 'proposed')
),
lines2 AS (
    SELECT 
        osm_id as id, 
        way as geom 
    FROM osm_polygon
    WHERE 
      (tags @> hstore('junction','roundabout') 
      OR tags @> hstore('highway','mini_roundabout'))
      AND NOT (tags ? 'proposed') 
)
SELECT
    'Roundabout' as label, 
    ST_AsText((ST_Dump(ST_LineMerge(ST_Multi(ST_Collect(geom))))).geom ) as geom
FROM lines
UNION
SELECT
    id::text as label, 
    ST_AsText(ST_ExteriorRing(geom))
FROM lines2
LIMIT 9999

/*

Other solution attempts:

A. Most promising but complex and not yet working: WITH RECURSIVE http://workshops.boundlessgeo.com/postgis-intro/advanced_geometry_construction.html
http://blog.cleverelephant.ca/2010/07/network-walking-in-postgis.html

B. This could be still an approach to follow as part of WITH RECURSIVE :

(
    SELECT *, COUNT(1) FROM lines l1 
    JOIN lines l2 ON ST_Intersects(l1.geom, l2.geom) WHERE l1.id <> l2.id
    GROUP BY l1.id, l1.geom
) as tmp
WHERE count = 2

Source: http://gis.stackexchange.com/questions/147316/merging-linestrings-in-qgis-using-postgresql-and-postgis


C. This query works and is from OpenMapTiles and assumes lines to join have common "name, highway, ref" - which is not the optimal condition here: 

SELECT
    -- NOTE: Make sure no MULII... are generated
    (ST_Dump(geometry)).geom AS geometry,
    -- NOTE: The osm_id is no longer the original one which can make it difficult
    -- to lookup road names by OSM ID
    member_osm_ids[0] AS osm_id,
    member_osm_ids,
    name,
    ref,
    highway,
    z_order
FROM (
    SELECT
        ST_LineMerge(ST_Collect(way)) AS geometry,
        name,
        ref,
        highway,
        min(z_order) AS z_order,
        array_agg(DISTINCT osm_id) AS member_osm_ids
    FROM osm_line
    -- NOTE: We only care about highways (not railways) for labeling
    WHERE (name <> '' OR ref <> '') AND NULLIF(highway, '') IS NOT NULL
    GROUP BY name, highway, ref
) AS highway_union

*/ 

Alle Geonamen, die bereits in Schweizerdeutscher Mundart ("name:gsw") erfasst sind

 SELECT 
   ST_AsText(way) geom,
   tags->'name:gsw' AS label
 FROM osm_poi
 WHERE tags ? 'name:gsw'

Finding a string part of a tag key or a tag value

 SELECT osm_id, name, key
 FROM (SELECT osm_id, name, skeys(tags) AS key 
       from osm_point) tmp
 WHERE key LIKE 'generator%';
 SELECT osm_id, name, value
 FROM (SELECT osm_id, name, svals(tags) AS value 
       FROM osm_polygon) tmp
 WHERE value LIKE 'nuclear%';

Alle Zoos der Schweiz

(Tipp: http://labs.geometa.info/postgisterminal/?xapi=*[tourism=zoo]):

  SELECT ST_AsText(way) AS geom, name||' '||osm_id AS label
  FROM osm_all
  WHERE tags @> hstore('tourism','zoo')  

Schweizer Kernkraftwerke mit 40 Km-Puffer

  SELECT
   ST_AsText(ST_Buffer(ST_Centroid(way),40000)) AS geom,
   name AS label
  FROM osm_all
  WHERE tags @> hstore('generator:source','nuclear')  

Alle Restaurants mit Namen 'Rössli' der Schweiz

  SELECT ST_AsText(way) AS geom, name AS label
  FROM osm_point
  WHERE amenity = 'restaurant'
  AND name ILIKE '%rössli%'  

Alle 4000er Berggipfel der Schweiz

  CREATE OR REPLACE FUNCTION as_numeric(text) RETURNS NUMERIC AS $$
  -- Inspired by http://stackoverflow.com/questions/16195986/isnumeric-with-postgresql/16206123#16206123
  DECLARE test NUMERIC;
  BEGIN
       test = $1::NUMERIC;
       RETURN test;
  EXCEPTION WHEN others THEN
       RETURN -1;
  END;
  $$ STRICT
  LANGUAGE plpgsql IMMUTABLE;

  SELECT ST_AsText(way) AS geom, osm_id||', '||coalesce(name,'No name')||', '||as_numeric(ele)::int||'m' AS label
  FROM osm_point
  WHERE tags @> 'natural=>peak' AND as_numeric(ele) between 4000 and 9999
  ORDER BY as_numeric(ele) desc;

  oder:

  SELECT ST_AsText(way) AS geom, name||','||ele AS label
  FROM osm_point
  WHERE "natural" = 'peak'
  AND  to_number('0'||ele, '99999999999.000')::int >= 4000
  
  oder:
  
  SELECT ST_AsText(way) AS geom, COALESCE(name, '')||' '||osm_id AS label
  FROM osm_poi
  WHERE tags @> hstore('natural', 'peak')
  AND CAST(regexp_replace(hstore("tags")->'ele', '[^0-9\.]', '', 'g') AS real) >= 4000  

Alle Seen der Schweiz

Gemäss https://twitter.com/MySwitzerland_e/status/923576854400897024 sind es über 1500. Nach internationalen Kriterien soll ein See mind. 1 Hectar (10000m2) sein. Der Wert scheint in der Schweiz etwas schmaler zu sein (z.B. Anenseeli 900m2?).

 SELECT ST_AsText(ST_PointOnSurface(way)) AS geom, name AS label
 FROM osm_polygon
 WHERE ( 
   tags @> hstore('water','lake') 
   OR (tags->'natural') IN ('water','lake')
 )
 AND ST_Area(way) > 200
 AND name IS NOT NULL

Alle Aussichtspunkte im Kanton Zürich, die höher als 500 m ü.M. sind

  SELECT ST_AsText(a.way) geom, COALESCE(name, '')||' '||ele||' m ü.M' AS label
  FROM 
    osm_point AS a, 
    (SELECT way FROM osm_polygon 
     WHERE name = 'Zürich'
     AND tags @> hstore('admin_level','4')) AS b
  WHERE ST_Contains(b.way,a.way)
  AND tags @> hstore('tourism','viewpoint')
  AND  to_number('0'||ele, '99999999999.000')::int >= 500   

Alle Picnic-Plätze und Aussichtspunkte im aktuellen Kartenausschnitt

  SELECT ST_AsText(way) AS geom, name AS label
  FROM osm_point
  WHERE tourism IN ('picnic_site','viewpoint')
  AND ST_Contains(mapextent(), way)  

Alle Schulhäuser im Umkreis von 40 km aller Kernkraftwerke

  SELECT ST_AsText(a.way) AS geom, '' AS label
  FROM 
    osm_poi AS a, 
    (SELECT ST_Buffer(ST_Centroid(way),40000) AS way
     FROM osm_poi
     WHERE tags @> hstore('generator:source','nuclear')) AS b
  WHERE ST_Within(a.way,b.way)
    AND a.tags @> hstore('amenity', 'school')
  UNION
  SELECT ST_AsText(ST_Buffer(ST_Centroid(way),40000)) AS geom, COALESCE(name, '') AS label
  FROM osm_poi
  WHERE tags @> hstore('generator:source','nuclear')
  LIMIT 2000  

Die nächsten 10 Bars in der Nähe von 'mylocation' (ungeachtet der Distanz)

 
  SELECT ST_AsText(osm_poi.way) AS geom, COALESCE(name, '') AS label
  FROM 
    osm_poi, 
    (SELECT ST_Transform(ST_GeomFromText('POINT(8.81638 47.22666)', 4326), 3857) AS way) AS mylocation
  WHERE osm_poi.tags @> hstore('amenity', 'bar')
  ORDER BY osm_poi.way <-> mylocation.way
  LIMIT 10  

or

 
  WITH mylocation AS ( 
    SELECT 
     ST_Transform(ST_GeomFromText('POINT(7.57369 47.54609)', 4326), 3857) AS way,
     'You are here'::text as label
  ),
  mypois AS (
    SELECT 
      osm_poi.way, 
      COALESCE(name, '*') AS label
    FROM osm_poi, mylocation
    WHERE osm_poi.tags @> hstore('historic', 'castle')
    ORDER BY osm_poi.way <-> mylocation.way
    LIMIT 10
  )
  SELECT 
    ST_AsText(way) AS geom, 
    label
  FROM mypois
  UNION ALL
  SELECT 
    ST_AsText(way) AS geom, 
    label
  FROM mylocation  

Die nächsten 100 Restaurants in der Nähe von 'mylocation' im Umkreis von max. 20 Km (Luftlinie)

  SELECT ST_AsText(osm_poi.way) AS geom, COALESCE(name, '') AS label
  FROM 
    osm_poi, 
    (SELECT ST_Transform(ST_GeomFromText('POINT(8.81638 47.22666)', 4326), 3857) AS way) AS mylocation
  WHERE ST_DWithin(osm_poi.way, mylocation.way, 20000)
  AND osm_poi.tags @> hstore('amenity', 'restaurant')
  ORDER BY osm_poi.way <-> mylocation.way
  LIMIT 100  

Alle Strassen mit Namen "Bahnhofstrasse" im Kt.ZH

 
  SELECT ST_AsText(ST_LineMerge(w.way))
  FROM osm_line w, 
       (SELECT way FROM osm_polygon 
        WHERE name = 'Zürich'
        AND tags @> hstore('admin_level','4')) AS p
  WHERE w.name = 'Bahnhofstrasse'
  AND (hstore("tags")->'highway') IS NOT NULL
  AND ST_Within(w.way, p.way)   

Bounding Box gegeben ein geografischer Name (Beispiel "Entlebuch")

  • Enter Entlebuch in http://boundingbox.klokantech.com/
  • Select DublinCore (and reduce to 4 digits) => westlimit=7.9891; southlimit=46.9419; eastlimit=8.3543; northlimit=47.1263
  • Place coords. in Polygon:
 
  SELECT ST_AsText( ST_Transform( ST_GeomFromText(
    'MULTIPOLYGON(((7.9891 47.1263, 8.3543 47.1263, 
                    8.3543 46.9419, 7.9891 46.9419, 
                    7.9891 47.1263)))', 4326), 3857) )   

Alle OeV-Haltestellen

  • Es hat 4000(!) mehr als in OSM als Didok/uic_name eingetragene Nodes?
  • Tags-Zusammenstellung aus DIDOK: https://github.com/datendelphin/didok/blob/master/db-import/import_stops.py
  • "SELECT way FROM osm_polygon WHERE osm_id=-51701" gibt die Schweizer Grenze zurück (EOSMDBOne enthält auch Teile des angrenzenden Auslands!).
  • Performance Note: The presence of a subquery is blocking the inlining of st_within, so putting the subquery into the from-clause instead.
 -- Just run this e.g. in PostGIS Terminal:
  SELECT count(*) 
  FROM 
    osm_point AS osm,
    (SELECT way FROM osm_polygon WHERE osm_id=-51701) AS ch
  WHERE ST_Within(osm.way, ch.way)
  AND NOT (tags ? 'uic_name')
  AND (tags @> 'railway=>station' OR tags @> 'railway=>halt' OR tags @> 'railway=>tram_stop' 
       OR tags @> 'highway=>bus_stop' OR tags @> 'amenity=>bus_station' 
       OR tags @> 'amenity=>ferry_terminal' OR tags @> 'railway=>funicular' OR tags @> 'aerialway=>station')
  -- 5291, query runtime: 1913ms.

  -- Use PostGIS Terminal to visualize the stops not in Didok (uic_name):
  SELECT ST_AsText(way) as geom, coalesce(name,'') as label
  FROM 
    osm_point AS osm,
    (SELECT way FROM osm_polygon WHERE osm_id=-51701) AS ch
  WHERE ST_Within(osm.way, ch.way)
  AND NOT (tags ? 'uic_name')
  AND (tags @> 'railway=>station' OR tags @> 'railway=>halt' OR tags @> 'railway=>tram_stop' 
       OR tags @> 'highway=>bus_stop' OR tags @> 'amenity=>bus_station' 
       OR tags @> 'amenity=>ferry_terminal' OR tags @> 'railway=>funicular' OR tags @> 'aerialway=>station')

  -- Get result as GeoCSV/WKT - to be run e.g. in psql or pgAdmin and export result to file:
  -- Round to mm precision (0.00000001 fractional part) from float using SnapToGrid. 
  -- See https://en.wikipedia.org/wiki/Decimal_degrees 
  FROM 
    osm_point AS osm,
    (SELECT way FROM osm_polygon WHERE osm_id=-51701) AS ch
  WHERE ST_Within(osm.way, ch.way)
  AND NOT (tags ? 'uic_name')
  AND (tags @> 'railway=>station' OR tags @> 'railway=>halt' OR tags @> 'railway=>tram_stop' 
       OR tags @> 'highway=>bus_stop' OR tags @> 'amenity=>bus_station' 
       OR tags @> 'amenity=>ferry_terminal' OR tags @> 'railway=>funicular' OR tags @> 'aerialway=>station')
  -- Resultat in https://gist.github.com/sfkeller/a18675d2190c4e121549e745474dc80f    

Die nächsen 15 Berge rund um Zürich

  • Next 15 mountains near Zurich
  • Kriterien: Höher als 1111m und mit Wikipedia-Tag (und natürlich mit Namen).
 SELECT 
    COALESCE(name, '') AS name,
    to_number('0'||ele, '99999999999.000')::int AS height,
    round(ST_Distance(start.geom, way))::int AS distance,
    CASE 
      WHEN round(degrees(ST_Azimuth(start.geom, way))) BETWEEN 315 AND 359 THEN 'north'
      WHEN round(degrees(ST_Azimuth(start.geom, way))) BETWEEN   0 AND  45 THEN 'north'
      WHEN round(degrees(ST_Azimuth(start.geom, way))) BETWEEN  45 AND 135 THEN 'east'
      WHEN round(degrees(ST_Azimuth(start.geom, way))) BETWEEN 135 AND 225 THEN 'south'
      WHEN round(degrees(ST_Azimuth(start.geom, way))) BETWEEN 225 AND 315 THEN 'west'
    END as direction,
    'https://de.wikipedia.org/wiki/'||(tags->'wikipedia') AS wikipedia,
    osm_id,
    ST_AsText(ST_SnapToGrid(ST_Transform(way,4326),0.00000001)) AS wkt 
  FROM 
    osm_point peak, 
    (SELECT ST_Transform(ST_GeomFromText('POINT(8.5393635 47.3781008)', 4326), 3857) AS geom) AS start
  WHERE tags @> hstore('natural', 'peak')
  AND to_number('0'||ele, '99999999999.000')::int >= 1111
  AND NULLIF(name, '') > ''
  AND tags ? 'wikipedia'
  ORDER BY way <-> start.geom
  LIMIT 15; 


Adresse eines Shops

  • Adresse eines Shops mit Attributen von seinen Umliegenden Objekten, wenn er selber keines hat.
  • Hier am Beispiel von Bankomaten (ATM):
  select 
    atm.osm_id, 
    atm.name, 
    atm.tags, 
    coalesce((atm.tags->'addr:street'), building.street, anyobject.street) as street, 
    coalesce((atm.tags->'addr:housenumber'), building.housenumber, anyobject.housenumber) as housenumber, 
    st_transform(atm.way,4326) as geom
  from osm_point as atm
  left join lateral (
    select (tags->'addr:street') as street, (tags->'addr:housenumber') as housenumber
    from osm_polygon
    where (tags->'building') is not null and (tags->'addr:street') is not null
    and st_within(atm.way, way)
  ) as building on true
  left join lateral (
    select (tags->'addr:street') as street, (tags->'addr:housenumber') as housenumber
    from osm_polygon
    where (tags->'addr:street') is not null
    order by atm.way <-> way
    limit 1
  ) as anyobject on true
  where tags @> 'amenity=>atm'::hstore
  order by street; 

Special SQL Queries

Extracting numbers out of OSM values

See also regex above...

 with tmp(txt) as (
   values ('1'),('2 '),('9.999999'),('12345678901.999'),('5.5'),('0.111')
     ,('-1.1'),('-0.00000001'),('-9999')
     ,(' 5.65 '),('xxx')
 )
 select 'Numeric', txt, coalesce(substring(txt, '([-]?0\.\d*[1-9]\d*|[-]?[1-9]\d*(\.\d+)?)')::numeric, 0)
 from tmp
 union
 select 'Integer', txt, coalesce(substring(txt, '([-]?[0-9]{1,19})')::bigint, 0)
 from tmp
 union
 select 'Number ', txt, to_number('0'||txt, 'S99999999999.000')::numeric
 from tmp
 order by 1,2

Regex: ....

Name string handling

Select names of all OSM-objects containing 'zoo' at start middle or end (using wildcard '%' in String):

 SELECT ST_AsText(way) geom,
 COALESCE(name, ' ')||' '||osm_id label
 FROM osm_all
 WHERE hstore(tags)->'name' ILIKE '%zoo%'

Select all parking lots for disabled persons (Rollstuhlparkplatz / wheel parking) within visible map area (without FIXME). Note: Tag capacity:disabled currently occurs always together with amenity={parking, parking_space, parking_entrance}):

 SELECT ST_AsText(way) geom, 
        COALESCE(name,' ')||' ('||COALESCE(SUBSTRING((tags->'capacity:disabled') FROM E'[0-9]+'),' ')||')' label
 FROM osm_poi
 WHERE (tags @> '"capacity:disabled"=>"yes"' 
        OR SUBSTRING((tags->'capacity:disabled') FROM E'[0-9]+')::int > 0)
 AND ST_Contains(mapextent(), way)

Select all motorways (higways) with speed limit greater or equal than 100:

 SELECT ST_AsText(way) geom
 FROM osm_line
 WHERE tags @> '"highway"=>"motorway"'
 AND COALESCE(SUBSTRING((tags->'maxspeed') FROM E'[0-9]+')::int,0) >= 100
 -- 1806!

Count tags except xxx...:

 SELECT array_upper(%# tags, 1) AS "#tags", count(*) AS "count"
 FROM osm_point
 WHERE array_upper(%# tags, 1) > 0
 AND NOT EXISTS (SELECT skeys FROM skeys(tags) WHERE skeys LIKE 'xxx:%')
 GROUP BY 1
 ORDER BY 1


Point Clustering

Update: See also PostGIS' new ST_ClusterWithin function.

Point Clustering with round/truncate to grid and centroid. Parameter 9000 defines the grid width (meter) and is dependent from the units of the coordinate reference systems. Fast - reduces 13358 Restaurants to 1163 - but with a "visual grid effect".

 SELECT
   ST_AsText( ST_Centroid(ST_Collect(position)) ) AS geom,
   CASE COUNT(position)
     WHEN 1 THEN COALESCE(max(name), '1') ELSE COUNT(position)::text
   END AS label,
   array_agg(osm_id) AS list_of_ids 
 FROM (
   SELECT way AS position, name, osm_id 
   FROM osm_poi
   WHERE tags @> hstore('amenity','restaurant')
 ) tmp
 GROUP BY ST_SnapToGrid(position, 9000);


Filter multivalued tags from OpenStreetMap

Transform a tag with multi-values (";" / semi-colon separated) - like 'cuisine' - to a clean, sorted text array! This allows for elegant filtering in WHERE clause. Array of text with a GIN index is the overall winner acording to this benchmark.

(Works with EOSMDBone, but not in the PostGIS Terminal, because of the ";" anti-db-injection function)

  with tmp as ( 
    select 
      osm_id, 
      name, 
      way, 
      lower(regexp_replace((tags->'cuisine'),'( )|;$','','g')) as cuisine  -- trim multi value string from spaces
    from osm_poi  -- view containing osm_point and osm_polygon with centroid.
    where (tags->'cuisine') is not null
  ),
  tmp2 as (
    select 
      osm_id,
      name, 
      way,  
      (select array(select unnest(regexp_split_to_array(cuisine,';')::text[]) order by 1)) as cuisine -- convert to sorted array
    from tmp
  )
  select cuisine, count(*), min(osm_id) from tmp2
  where cuisine @> '{swiss}' -- does cuisine-array contain 'swiss'? 
  --where cuisine @> '{swiss,regional}' -- matches "swiss AND regional" but also "regional AND swiss"
  group by cuisine
  order by 2 desc;

  /* Data Output:
  {swiss}                              51 33117844
  {regional,swiss}                      2 400694956
  {italian,swiss}                       2 354543026
  {breakfast,coffee,regional,swiss}     1 99418980
  {breakfast,regional,swiss}            1 2371979274
  {fondue,raclette,swiss}               1 380852254
  {italian,pasta,pizza,swiss}           1 5186344683
  {french,swiss}                        1 988716374
  {european,mediterranean,swiss,tapas}	1 3145378112
  */

This is some test code showing step-by-step how the multivalue string is being converted to an array:

 with osm_poi_raw (id, name, cuisine_str) as ( values 
    (1, 'Reschti 1', 'italian;swiss;pizzeria'), 
    (2, 'Reschti 2', 'french;swiss'),
    (3, 'Reschti 3', 'fast_food;american'),
    (4, 'Test 1', 'ok_1;ok_2; not_1;not_2 ; not_3 ;'),
    (5, 'Test 2', 'ok_1;ok_2;_not_1;not_2_;_not_3_;')
  ),
  osm_poi (id, name, cuisine_arr) as (
    --select id, name, regexp_split_to_array(cuisine_str,';') 
    --select id, name, regexp_split_to_array(regexp_replace(cuisine_str,'( |_)|;$','','g'),';') 
    select id, name, (select array(select unnest(regexp_split_to_array(regexp_replace(cuisine_str,'( |_)|;$','','g'),';')) order by 1))
    from osm_poi_raw
  )
  select * from osm_poi
  where (cuisine_arr @> '{swiss,french}' or cuisine_arr @> '{american}');

Marker Queries

marker query (must have exactly field names lon, lat, title and description):

 SELECT X(p2.way) AS lon, Y(p2.way) AS lat, 'Briefkasten' AS title, p2.ref AS description
 FROM planet_osm_polygon p1
 JOIN planet_osm_point p2 ON CONTAINS(p1.way, p2.way)
 WHERE p1.name = 'Uster'
 AND p2.amenity = 'post_box'

Extension with user defined marker icon (attribute 'icon'):

 SELECT X(p2.way) AS lon, Y(p2.way) AS lat, 'Briefkasten' AS title, p2.ref AS description, 'http://myserver/marker.png' as icon 
 FROM planet_osm_polygon p1
 JOIN planet_osm_point p2 ON CONTAINS(p1.way, p2.way)
 WHERE p1.name = 'Uster'
 AND p2.amenity = 'post_box'

Mouse Queries

Finding Nearest Point to Linestrings (Snapping)

  • Lösung mit Distanzen innerhalb 50m (siehe ST_DWithin).
  • Anwendung: PostGIS-Terminal aufrufen, Copy&Paste dieser Query, dann in Karte klicken.
  • Bemerkungen:
    • (_mouse_x, _mouse_y) liefern die Mauskoordinaten (anstelle mypos).
    • Die Lösung verwendet WITH-Klausel, da (_mouse_x, _mouse_y) vom Terminal nur einmal geparst werden (leider).
    • (tags->'highway') liefert die Strassenart, z.B. unclassified, pedestrian, etc.
    • (SELECT geom FROM tmp1) liefert die Mauskoordinaten
    • COALESCE(name,'NN') liefert NN wenn kein Name vorhanden
 WITH tmp1 AS ( 
   SELECT ST_SetSRID(ST_MakePoint(_mouse_x, _mouse_y), 3857) AS geom
 ), 
 tmp2 AS ( 
   SELECT osm_id, way, name, tags, ST_Distance((SELECT min(geom)::geometry FROM tmp1), way) as distance
   FROM osm_line
   WHERE (tags->'highway') IS NOT NULL
   AND ST_DWithin((SELECT geom FROM tmp1), way, 50)
 )
 SELECT 
   ST_AsText(way) AS geom, 
   COALESCE(name,'NN')||' ('||(tags->'highway')||')' AS label
 FROM tmp2
 ORDER BY distance
 LIMIT 1

(Fast) dasselbe wie oben aber mit KNN Index (<->) nun aber mit den 10 nächsten Linien (k-Nearest Neighbor (KNN) Index)

  • geordnet nach dem Zentrum der Bounding Boxes ('<->' Operator).
  • Hinweis: der <#> ordnet nach den Bounding Boxen selber (ab PostGIS 2.0 und PostgreSQL 9.1)
 WITH tmp1 AS ( 
   SELECT ST_SetSRID(ST_MakePoint(_mouse_x, _mouse_y), 3857) AS geom
 ), 
 tmp2 AS ( 
   SELECT osm_id, way, name, tags, ST_Distance((SELECT min(geom)::geometry FROM tmp1), way) as distance
   FROM osm_line
   WHERE (tags->'highway') IS NOT NULL
   ORDER BY way <-> (SELECT min(geom)::geometry FROM tmp1)
   LIMIT 10
 )
 SELECT 
   ST_AsText(way) AS geom, 
   COALESCE(name,'NN')||' ('||(tags->'highway')||')' AS label
 FROM tmp2
 ORDER BY distance
 LIMIT 1

Statistics

Youngest node in the database

 select 
   osm_id, 
   name, 
   (tags->'osm_timestamp')::timestamptz as timestamp, 
   (tags->'osm_version')::int as version,
   (tags->'osm_uid')::int as uid, 
   (tags->'osm_changeset')::int as changeset, 
   way
 from osm_point
 where osm_id=(select max(osm_id) from osm_point)


Keys containing problematic colons

List top 20 keys which contain problematic colons (":") and have value "yes", grouped by occurrence:

 WITH tmp AS (
   SELECT 
     SUBSTRING(key FROM '#"%#":%' for '#') AS key,
     value
   FROM ( 
     SELECT 
       osm_id, 
       ((each(tags)).key) AS key,
       ((each(tags)).value) AS value
     FROM osm_point
   ) AS osm_point
   WHERE value IN ('yes','other')
   AND key LIKE '%:%' 
 )
 SELECT key, count(*) AS count
 FROM tmp
 GROUP BY key
 ORDER BY count DESC, key ASC
 LIMIT 20

All tupels in all tables

 SELECT '1. '||to_char(count(*), '999G999G999')||' osm_point(s)'
 FROM osm_point
 UNION
 SELECT '2. '||to_char(count(*), '999G999G999')||' osm_poi(s)'
 FROM osm_poi
 UNION
 SELECT '3. '||to_char(count(*), '999G999G999')||' osm_line(s)'
 FROM osm_line
 UNION
 SELECT '4. '||to_char(count(*), '999G999G999')||' osm_polygon(s)'
 FROM osm_polygon
 UNION
 SELECT '5. '||to_char(count(*), '999G999G999')||' osm_nodes'
 FROM osm_nodes
 UNION
 SELECT '6. '||to_char(count(*), '999G999G999')||' osm_ways'
 FROM osm_ways
 UNION
 SELECT '7. '||to_char(count(*), '999G999G999')||' osm_rels'
 FROM osm_rels
 ORDER BY 1

All Tag-Value-Pairs

Kann anstelle mit osm_point auch mit osm_all durchgeführt werden.

 --
 -- Key-Value Statistics of OSM Data
 -- Return all key-value-pairs of type 'enum' without some 
 -- types numeric, date/time etc. chosen by hand:
 -- 
 -- Alternative with separate tag-value: 
 -- SELECT tmp.key as tag, tmp.value as value, count(*)::integer as freq
 SELECT tmp.key||'='||tmp.value as kvp, count(*)::integer as freq 
 FROM (
   SELECT (each(tags)).key as key, (each(tags)).value as value 
   FROM osm_point) as tmp
 WHERE (trim(value) !~ '^[-]*[0-9,.:\ ]+[m]*$')
 AND NOT (value ILIKE '%fixme%' OR key ILIKE '%fixme%')
 AND key NOT LIKE '%:%'
 AND key NOT LIKE '%description%'
 AND key NOT LIKE '%comment%'
 AND key NOT LIKE '%name'
 AND key NOT LIKE 'uic_%'
 AND key NOT LIKE '%ref'
 AND key NOT ILIKE '%fixme%'
 AND key NOT ILIKE '%todo%'
 AND key NOT IN ('osm_timestamp', 'name','operator','_picture_','_waypoint_','address','alt','is_in','url','website','wikipedia','email',
                 'converted_by','phone','information','opening_hours','date','time','collection_times','colour','fee',
                 'population','access','noexit','towards','bus_routes','busline','lines','type','denotation',
                 'CONTINUE','continue','copyright','stop')
 GROUP BY tmp.key, tmp.value
 HAVING COUNT(*) > 1
 ORDER by freq DESC; -- ca. 500 sec.!