09. PostGIS - prostorové SQL

Přednáška

Úvod do PostGIS

-- Geometrické konstruktory

SELECT ST_GeomFromText('POINT(-686651 -1058147)');

SELECT ST_Point(-686651, -1058147);

SELECT 'POINT(-686651 -1058147)'::geometry;

SELECT ST_GeomFromText('POINT(-686651 -1058147)', 5514);

SELECT ST_SetSRID(ST_Point(-686651, -1058147), 5514);

SELECT ST_GeomFromEWKT('SRID=5514;POINT(-686651 -1058147)');

SELECT 'SRID=5514;POINT(-686651 -1058147)'::geometry;

-- Prostorový index

SELECT count(*)
FROM   ruian.obce AS obce
JOIN   osm.zeleznice AS zelez
ON     ST_Intersects(zelez.geom, obce.geom);

-- Query planner
EXPLAIN
SELECT  count(*)
FROM    ruian.obce AS obce
JOIN    osm.zeleznice AS zelez
ON      ST_Intersects(zelez.geom, obce.geom);

-- ověření
SELECT   obce.kod, count(*) as pocet
FROM     ruian.obce AS obce
JOIN     osm.zeleznice AS zelez
ON       ST_Intersects(zelez.geom, obce.geom)
GROUP BY obce.kod
ORDER BY pocet DESC
LIMIT    10;

-- možná řešení
SELECT  count(DISTINCT obce.kod)
FROM    ruian.obce AS obce
JOIN    osm.zeleznice AS zelez
ON      ST_Intersects(zelez.geom, obce.geom);

-- explain analyze
EXPLAIN ANALYZE 
SELECT  count(DISTINCT obce.kod)
FROM    ruian.obce AS obce
JOIN    osm.zeleznice AS zelez
ON      ST_Intersects(zelez.geom, obce.geom);

EXPLAIN ANALYZE 
SELECT  count(DISTINCT obce.kod)
FROM    ruian.obce AS obce
JOIN    osm.zeleznice AS zelez
ON      zelez.geom && obce.geom
AND     _st_intersects(zelez.geom, obce.geom);

-- bez prostoroveho indexu
EXPLAIN ANALYZE 
SELECT  count(DISTINCT obce.kod)
FROM    ruian.obce AS obce
JOIN    osm.zeleznice AS zelez
ON     _st_intersects(zelez.geom, obce.geom);

Cvičení

Data Manupulation Language

Prostorové operátory

Vyhledání prvků v okolí daného bodu.

SELECT osm_id, st_astext(geom)
FROM   osm.pozarni_stanice
WHERE  geom @ st_expand(
 st_geomfromtext('POINT(-768581 -1063422)', 5514),
 20000);

-- vs.

SELECT osm_id, st_astext(geom)
FROM   osm.pozarni_stanice
WHERE  geom @ st_buffer(
 st_geomfromtext('POINT(-768581 -1063422)', 5514),
 20000);

-- vs.

SELECT osm_id, st_astext(geom)
FROM   osm.pozarni_stanice
WHERE  st_within(geom, st_buffer(
 st_geomfromtext('POINT(-768581 -1063422)', 5514),
 20000));

-- vs.

SELECT osm_id, st_astext(geom)
FROM   osm.pozarni_stanice
WHERE  geom <#> 'SRID=5514;POINT(-768581 -1063422)'::geometry < 20000;

Vyhledání prvků ležících kompletně/částečně uvnitř MMO.

SELECT   count(*)
FROM     osm.silnice AS s
JOIN
(
SELECT   geom
FROM     ruian.vusc
WHERE    nutslau = 'CZ042'
) AS u
ON       s.geom @ u.geom;

-- vs.

SELECT   count(*)
FROM     osm.silnice AS s
JOIN     ruian.vusc as u
ON       u.nutslau = 'CZ042'
AND      s.geom && u.geom;

Prostorová analýza

Vyhledání nejjižněji položeného zájmového bodu.

SELECT   osm_id, st_y(geom) AS y
FROM     osm.pozarni_stanice
ORDER BY y ASC
LIMIT    1;

Plošný prvek s vnitřní hranicí.

SELECT   nazev, st_area(geom)/1e6 AS area_km
FROM     ruian.vusc
WHERE    st_nrings(geom) > 1
ORDER BY area_km DESC LIMIT 1;

Centroid (referenční bod) plošného prvku.

SELECT nazev, st_astext(st_pointonsurface(geom))
FROM   ruian.vusc;

Určení výměry plošného prvku.

SELECT nazev, round(st_area(geom)/1e6) AS area_km
FROM   ruian.vusc;

Seřazení plošných prvků podle výměry.

SELECT   nazev, st_area(geom) AS plocha
FROM     ruian.vusc
ORDER BY plocha DESC;

Celková délka liniových prvků.

SELECT   round(sum(st_length(geom))/1000) AS dalnice_km
FROM     osm.silnice 
WHERE    typ = 1;

Plošná zakulacenost prvků (poměr kvadrátu obvodu vůči ploše).

SELECT nazev, (st_perimeter(geom) * st_perimeter(geom))
       / st_area(geom) AS hodnota
FROM    ruian.vusc;

Najít zájmové body, které se nacházejí do vzdálenosti od plošného prvku.

SELECT osm_id FROM osm.pozarni_stanice AS p
JOIN   ruian.obce AS k
ON     k.nazev = 'Kladno'
AND    st_distance(p.geom, k.geom) < 10000;

Zájmové body, které se nacházejí ve vzdálenosti od nejdelšího úseku liniového prvku.

SELECT   osm_id, st_astext(geom)
FROM     osm.pozarni_stanice
WHERE    st_distance(
(
SELECT   geom
FROM     osm.silnice
WHERE    typ = 1 
ORDER BY st_length(geom) DESC
LIMIT    1
),       geom) < 10000;

-- vs.

WITH s AS (
SELECT   geom
FROM     osm.silnice
WHERE    typ = 1 
ORDER BY st_length(geom) DESC
LIMIT    1)
SELECT   osm_id, st_astext(p.geom)
FROM     osm.pozarni_stanice as p
JOIN     s
ON       s.geom <-> p.geom < 10000;

Určit vzdálenost mezi zájmovými bodovými prvky.

SELECT (st_distance(
(
SELECT   geom
FROM     osm.pozarni_stanice
ORDER BY st_x(geom) DESC
LIMIT 1
), 
(
SELECT geom
FROM    osm.pozarni_stanice
ORDER BY st_x(geom) ASC
LIMIT 1
))/1e3
)::int as vzdalenost_km;

Prostorové predikáty

Vypsání sumárních údajů.

SELECT   k.kod, k.nazev,
 SUM(st_length(z.geom))/1000 AS zel_km
FROM     osm.zeleznice AS z
JOIN     ruian.vusc AS k
ON       st_contains(k.geom, z.geom) 
GROUP BY k.kod, k.nazev
ORDER BY zel_km;

-- vs.

SELECT   k.kod, k.nazev,
 SUM(st_length(z.geom))/1000 AS zel_km
FROM     osm.zeleznice AS z
JOIN     ruian.vusc AS k
ON       st_intersects(z.geom, k.geom) 
GROUP BY k.kod, k.nazev
ORDER BY zel_km;

-- vs.

SELECT   k.kod, k.nazev,
 SUM(st_length(st_intersection(z.geom, k.geom)))/1000 AS zel_km
FROM     osm.zeleznice AS z
JOIN     ruian.vusc AS k
ON       z.geom && k.geom
GROUP BY k.kod, k.nazev
ORDER BY zel_km;

Liniové prvky, které leží do vzdálenosti od zájmových bodových prvků.

SELECT count(DISTINCT s.osm_id)
FROM   osm.silnice AS s
JOIN   osm.pozarni_stanice AS p
ON     st_dwithin(s.geom, p.geom, 300);

-- vs.

SELECT count(DISTINCT s.osm_id)
FROM   osm.silnice AS s
JOIN   osm.pozarni_stanice AS p
ON     p.geom && st_expand(s.geom, 300)
AND    p.geom <-> s.geom < 300;

-- vs.

WITH pb AS (
SELECT st_buffer(geom, 300, 120) as geom
FROM   osm.pozarni_stanice
)
SELECT count(DISTINCT s.osm_id)
FROM   osm.silnice AS s
JOIN   pb
ON     st_intersects(s.geom, pb.geom);

Vytvořit liniové prvky, které leží do vzdálenosti od zájmových bodových prvků.

WITH pb AS (
SELECT st_buffer(geom, 300) as geom
FROM   osm.pozarni_stanice
)
SELECT DISTINCT s.osm_id,
 st_length(st_intersection(s.geom, pb.geom)) as geom
FROM   osm.silnice AS s
JOIN   pb
ON     st_intersects(s.geom, pb.geom) ORDER BY osm_id LIMIT 5;

-- vs.

WITH pb AS (
SELECT st_buffer(geom, 300, 42) as geom
FROM   osm.pozarni_stanice
)
SELECT DISTINCT s.osm_id,
 st_length(st_intersection(s.geom, pb.geom)) as geom
FROM   osm.silnice AS s
JOIN   pb
ON     st_intersects(s.geom, pb.geom) ORDER BY osm_id LIMIT 5;

Prostorová agregace

Najdi průsečíky železnic, kde se protínají přesně čtyři linie.

SELECT ST_AsText(point) FROM
 (SELECT ST_StartPoint(geom) AS point FROM osm.zeleznice
  UNION ALL
  SELECT ST_EndPoint(geom) AS point FROM osm.zeleznice)
AS a GROUP BY point HAVING COUNT(*) = 4;

Konverze prvků

CREATE TABLE obce_s AS
 SELECT kod, nazev, ST_GeometryN(geom,
  generate_series(1, ST_NumGeometries(geom))) AS geom
  FROM ruian.obce;

-- vs.

 CREATE TABLE obce_s AS
 SELECT kod, nazev, (ST_Dump(geom)).geom as geom
  FROM ruian.obce;

Konverze prvku (singleparts) na multiprvek (multiparts).

CREATE TABLE obce_m AS
 SELECT   kod, nazev, ST_Union(geom) AS geom
 FROM     obce_s
 GROUP BY kod, nazev;

Rozdělení linie na daný počet segmentů.

WITH linie AS (
    SELECT st_boundary(geom) AS geom
    FROM ruian.staty
), a AS (
   SELECT generate_series(0,999,1) AS a
)
SELECT ST_LineMerge(
 ST_LineSubstring(
  geom, a * 0.001, (a + 1) * 0.001)
) AS geom
FROM a, linie;

-- jako procedura...

CREATE OR REPLACE FUNCTION rozdel_pocet(
p_geometry GEOMETRY, p_pocet_segmentu INTEGER)
    RETURNS TABLE(v_hranice GEOMETRY)
AS $$
BEGIN
    RETURN QUERY
        WITH series AS (
            SELECT generate_series(
               0, p_pocet_segmentu-1, 1) AS series
        )
    SELECT ST_LineMerge(ST_LineSubstring(
        p_geometry, series * 0.001, (series + 1) * 0.001)
    ) AS v_hranice
    FROM series;
END; $$
LANGUAGE 'plpgsql';

SELECT rozdel_pocet(
 p_geometry => st_boundary(geom), p_pocet_segmentu => 1000
)
FROM ruian.staty;

Rozdělení linie na segmenty s danou délkou.

WITH linie AS (
   SELECT ogc_fid, st_boundary(geom) AS geom
   FROM ruian.staty
)
SELECT ST_LineSubstring(geom, 2000 * n / length,
  CASE
  WHEN 2000 * (n + 1) < length
  THEN 2000 * (n + 1) / length
  ELSE 1
  END) AS geom
FROM 
(SELECT ogc_fid,
  st_linemerge(geom) AS geom,
  st_length(geom) AS length
  FROM linie
  ) AS t
CROSS JOIN generate_series(0, 2000) AS n
WHERE n * 2000 / length < 1;

Data Definition Language

Vytvoření databáze.

CREATE DATABASE jmeno_db;
CREATE EXTENSION postgis;

Vytvoření tabulky (id, popisná data, geometrie).

CREATE TABLE poi (fid serial primary key,
                  nazev varchar(100),
                  geom geometry(point, 5514));

INSERT INTO poi (nazev, geom) VALUES (
 'FSv CVUT v Praze',
 st_transform(
  'SRID=4326;POINT (14.3880675 50.1040214)'::geometry,
 5514
));

SELECT * FROM geometry_columns WHERE f_table_name='poi';

Vytvoření tabulky jako výsledek dotazu:

CREATE TABLE okresy AS
 SELECT okreskod, st_union(geom) AS geom FROM ruian.obce
 GROUP BY okreskod;

-- Přidání primárního klíče:
ALTER TABLE okresy ADD COLUMN fid serial;
ALTER TABLE okresy ADD PRIMARY KEY (fid);

-- Omezení geometrie:
ALTER TABLE okresy ALTER COLUMN geom
  TYPE geometry(MULTIPOLYGON, 5514)
  USING st_multi(st_setsrid(geom, 5514));

--Sestavení prostorového indexu:
CREATE INDEX ON okresy USING gist (geom);

-- vs.
CREATE TABLE okresy (
fid serial PRIMARY KEY, okreskod int,
geom geometry(multipolygon, 5514)
);

INSERT INTO okresy (okreskod, geom)
 SELECT okreskod, st_multi(st_union(geom)) AS geom
 FROM ruian.obce GROUP BY okreskod;

CREATE INDEX ON okresy USING gist (geom);

Odstranění atributu geometrie/tabulky.

ALTER TABLE okresy DROP COLUMN geom;
DROP TABLE okresy;

-- vs.
SELECT DropGeometryColumn('okresy', 'geom');
SELECT DropGeometryTable('okresy');

Datová integrita

-- Validní prvek:
SELECT ST_IsValid('LINESTRING(0 0, 1 1)');

-- Invalidní prvek:
SELECT ST_IsValid('LINESTRING(0 0, 0 0)');
SELECT ST_IsValidReason('LINESTRING(0 0, 0 0)');

SELECT ST_IsSimple('LINESTRING(0 0, 1 1, 1 0, 0 1)');

SELECT ST_IsClosed('LINESTRING(0 0, 1 1, 1 0, 0 1, 0 0)');

SELECT ST_IsRing('LINESTRING(0 0, 1 0, 0 1, 0 0)');