PostGIS

Z FreeGIS portál
Přejít na: navigace, hledání
Postgis.png

PostGIS je rozšíření objektově-relačního databázového systému PostgreSQL pro podporu geografických objektů. PostGIS implementuje specifikaci Simple Features konsorcia Open Geospatial Consortium.

Rozšíření

Pro experimenty používejte cvičnou databázi.

Návod na instalaci PostGIS, poznámky k administraci.

Příklady

GIS 1

Úvod do zpracování prostorových dat

Vytvoření databáze

Pro verzi PostgreSQL 9.1 a vyšší

createdb <databáze>
psql -d <databáze> -c "CREATE EXTENSION postgis;"

Pro verzi PostgreSQL menší než 9.1

createdb <databáze>
createlang plpgsql <databáze>
psql -d <databáze> -f cesta/k/postgis.sql
psql -d <databáze> -f cesta/k/spatial_ref_sys.sql 

Poznámka: Pro nahrání dávky postgis.sql jsou vyžadována práva "superuživatele" (superuser) databázového systému PostgreSQL.

Import dat

V případě, že jsou atributová v jiném kódování, je nutné nastavit před importem proměnnou prostředí PGCLIENTENCODING, např.

export PGCLIENTENCODING=WIN1250

SQL

INSERT INTO zb (cat, zb_geom) VALUES (1, ST_GeomFromText('POINT(13.30150173 49.79076912)', 4326));

nebo

INSERT INTO zb (cat, zb_geom) VALUES (1, ST_SetSRID(ST_MakePoint(13.30150173,49.79076912), 4326));

Viz geometrické konstruktory.

ESRI Shapefile (shp2pgsql, ogr2ogr)

shp2pgsql -s 4326 -D -I -W latin2 cr.shp cr | psql pgis_yfsg
OGR včetně transformace S-JTSK → WGS-84
ogr2ogr -f PostgreSQL -s_srs '+proj=krovak \
+a=6377397.155 +rf=299.1528128 +no_defs \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +to_meter=1.0' \
-t_srs EPSG:4326 \
pg:dbname=pgis_yfsg cr.shp

GRASS GIS

Modul v.in.ogr.

v.in.ogr dsn="pg:dbname=pgis_yfsg" layer=obce format=PostgreSQL

Příklad skriptu pro import vektorových vrstev z databáze 'pgis_student' schéma 'gis1'.

#!/usr/bin/python                                                                                                                                                           

import sys
import psycopg2
import grass.script as grass

dbname = 'pgis_student'
schema = 'gis1'

def import_maps(maps):
    global dbname
    global schema
    for name in maps:
        grass.message("Importuji vektorovou mapu <%s>..." % name)
        grass.run_command('v.in.ogr',
                          quiet = True,
                          overwrite = True,
                          dsn = 'PG:dbname=%s' % dbname,
                          layer = schema + '.' + name,
                          output = name)
        grass.run_command('v.info',
                          flags = 't',
                          map = name)
def main():
    global dbname
    global schema
    conn = psycopg2.connect("dbname=%s" % dbname)
    cur = conn.cursor()

    cur.execute("SELECT tablename FROM pg_tables WHERE schemaname = '%s';" % schema)
    maps = list()
    for row in cur.fetchall():
        maps.append(row[0])

    cur.close()
    conn.close()

    import_maps(maps)

    return 0

if __name__ == "__main__":
    sys.exit(main())

Další formáty podporované knihovnou OGR (ogr2ogr)

GPX
ogr2ogr -f PostgreSQL -t_srs EPSG:4326 pg:dbname=pgis_yfsg yfsg-w.gpx waypoints

Poznámky

Rozsah dat

SELECT ST_Extent(the_geom) from points;
-- EPSG:4326 --
SELECT ST_Extent(ST_Transform(the_geom, 4326)) from points;

Self-intersection

Příklad detekce průniku polygonů lesních porostů (data OpenStreetMap).

CREATE VIEW lesy AS SELECT osm_id,way FROM czech_polygon WHERE landuse = 'forest';

CREATE AGGREGATE array_accum (
    sfunc = array_append,
    basetype = anyelement,
    stype = anyarray,
    initcond = '{}'
);

CREATE TABLE ilesy AS SELECT COUNT(DISTINCT g1.osm_id) AS count1, COUNT(DISTINCT g2.osm_id) AS count2,
        array_accum(DISTINCT g1.osm_id) AS array_accum1, array_accum(DISTINCT g2.osm_id) AS array_accum2,
        st_collect(DISTINCT g1.way) AS st_collect1,
        st_collect(DISTINCT g2.way) AS st_collect2,
        st_intersection(g1.way, g2.way), st_area(st_intersection(g1.way, g2.way))
FROM lesy g1, lesy g2
WHERE g1.osm_id < g2.osm_id AND st_intersects(g1.way, g2.way) AND
        st_isvalid(g1.way) AND st_isvalid(g2.way)
GROUP BY st_intersection(g1.way, g2.way)
ORDER BY COUNT(DISTINCT g1.osm_id), st_area(st_intersection(g1.way, g2.way));

Round()

Jaká je výměra lesa na mapovém listu 'M-33-63-A' v km2?

SELECT ROUND(SUM(ST_Area(ST_Intersection(kltm50.geom, lesy.geom)))/1e6)
FROM   kltm50
JOIN   lesy 
ON     kltm50.nazev = 'M-33-63-A'
AND    ST_Intersects(kltm50.geom, lesy.geom);
125
SELECT ROUND(CAST (SUM(ST_Area(ST_Intersection(kltm50.geom, lesy.geom)))/1e6) AS NUMERIC, 2)
FROM   kltm50
JOIN   lesy 
ON     kltm50.nazev = 'M-33-63-A'
AND    ST_Intersects(kltm50.geom, lesy.geom);
124.71

Prvek → multi-prvek

Příklad pro převod 'point' na 'multipoint' na základě kategorie (cat). Vstupní vrstva

SELECT cat,st_astext(wkb_geometry) from points;
 cat |                st_astext                 
-----+------------------------------------------
   1 | POINT(667291.931935888 278497.259214522)
   1 | POINT(387537.265811116 85478.9623985263)
   1 | POINT(667082.230368966 25935.0171310276)
   1 | POINT(272516.273379605 14288.4512718142)
   1 | POINT(469037.995578104 22458.0219971467)
   1 | POINT(299057.596718834 48860.4834410143)
   1 | POINT(705577.849591098 55690.8445368252)
   2 | POINT(357105.880943444 282025.922685119)
   2 | POINT(379192.913151337 240472.825207718)
   2 | POINT(512884.814365169 184280.16415778)
(10 rows)
SELECT cat,st_astext(st_union(wkb_geometry)) from points group by cat;
 cat | st_numgeometries 
-----+------------------
   1 |                7
   2 |                3
(2 rows)

Příklad pro tabulku obcí ze cvičné databáze PostGIS.

SET search_path TO public,gis1;

Počet obcí v ČR

SELECT count(DISTINCT kodob) from obce;
 count 
-------
  6249

Řada obcí je tvořena více polygony

SELECT nazev,count(*) from obce GROUP BY kodob,nazev HAVING count(*) > 1;
         nazev         | count 
-----------------------+-------
 Šternberk             |     2
 Mišovice              |     2
 Skuteč                |     2
 Krásensko             |     2
 Telč                  |     2
 Hořepník              |     2
 Nýrsko                |     2
 Přelouč               |     2
 Chlumec               |     2
 Spořice               |     2
 Kryry                 |     2
 Senice na Hané        |     2
 Vysoký Újezd          |     2
 Odry                  |     2
...
 Řehlovice             |     2
 Bohdalice-Pavlovice   |     2
 Letiny                |     2
 Bruntál               |     3
(95 rows)

Cílem je vytvořit novou tabulku obcí, kde pro tyto záznamy bude místo 'polygonu' použit 'multipolygon', tj. bude platit "jeden záznam == jedna obec".

-- create multipolygon layer
CREATE TABLE obce_tmp1 AS SELECT kodob,st_union(geom) AS geom FROM obce GROUP BY kodob;
-- copy original data
CREATE TABLE obce_tmp2 AS SELECT * FROM obce;
-- remove duplicate rows
DELETE FROM obce_tmp2 WHERE ctid NOT IN (SELECT MAX(ctid) FROM obce GROUP BY kodob);
-- replace polygon by multipolygon - geometry column
SELECT DropGeometryColumn('obce_tmp2', 'geom');
SELECT AddGeometryColumn('obce_tmp2', 'geom', 2065, 'MultiPolygon', 2);
-- copy geometry data (multipolygons)
UPDATE obce_tmp2 AS a SET geom = st_multi(b.geom) FROM obce_tmp1 AS b WHERE a.kodob = b.kodob;
-- clean up
DROP TABLE obce_tmp1;
ALTER TABLE obce_tmp2 RENAME TO obce_multi;
-- define primary key (required by QGIS)
ALTER TABLE obce_multi ADD primary key(ogc_fid);
-- build spatial index
CREATE INDEX obce_multi_geom_idx ON obce_multi using gist (geom);
-- update geom-related columns
UPDATE obce SET area = st_area(geom);
UPDATE obce SET perimeter = st_perimeter(geom);

Segmentace lomené čáry

create table line_test as select ST_GeomFromText('LINESTRING(0 0, 1 1, 2 2, 3 3, 4 4)') as geom;
  • podle lomových bodů
SELECT ST_AsText( ST_MakeLine(sp,ep) ) FROM
(
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 line_test
) AS foo;
 LINESTRING(0 0,1 1)
 LINESTRING(1 1,2 2)
 LINESTRING(2 2,3 3)
 LINESTRING(3 3,4 4)
  • podle vzdálenosti
alter table line_test add column gid serial;
CREATE OR REPLACE FUNCTION line_split(geometry, integer)
 RETURNS SETOF geometry AS 
$$ 
SELECT ST_Line_Substring(the_geom, $2 * n / length,
  CASE
	WHEN $2 * (n + 1) < length THEN $2 * (n + 1) / length
	ELSE 1
  END) AS the_geom
FROM 
(SELECT ST_LineMerge($1) AS the_geom,
  ST_Length($1) AS length) AS t
CROSS JOIN generate_series(0, $2) AS n
WHERE n * $2 / length < 1;
$$
LANGUAGE 'sql';
SELECT ST_AsText(line_split(geom, 1)) FROM line_test;

Alternativně:

select ST_AsText(ST_segmentize(geom, 1)) FROM line_test;
LINESTRING(0 0,
           0.707106781186548 0.707106781186548,
           1 1,
           1.70710678118655 1.70710678118655,
           2 2,
           2.70710678118655 2.70710678118655,
           3 3,
           3.70710678118655 3.70710678118655,
           4 4)
create view seg_line_test as select ST_segmentize(geom, 1) as geom FROM line_test;
SELECT ST_AsText( ST_MakeLine(sp,ep) ) FROM
(
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 seg_line_test
) AS foo;
 LINESTRING(0 0,0.707106781186548 0.707106781186548)
 LINESTRING(0.707106781186548 0.707106781186548,1 1)
 LINESTRING(1 1,1.70710678118655 1.70710678118655)
 LINESTRING(1.70710678118655 1.70710678118655,2 2)
 LINESTRING(2 2,2.70710678118655 2.70710678118655)
 LINESTRING(2.70710678118655 2.70710678118655,3 3)
 LINESTRING(3 3,3.70710678118655 3.70710678118655)
 LINESTRING(3.70710678118655 3.70710678118655,4 4)

Překrývající se polygony

Testovací data:

CREATE TABLE ref (
    cat integer,
    geom geometry(Polygon)
);

INSERT INTO ref VALUES ('010300002032BF0D0001000000110000001E59C611F3491D410A11BAB453351141EC499BDF68082041D8A65F4924731141348442B4C13521418B77E658D2FB0C41E8E43A937A9221419BF8A3719E080841348442B4C13521412F5D3BE506E50441A5FC2FF222102041AE669D717D4A0441E04A333FC1311E417EE97EED60BD05419D62289E69AF1C41C89B4A27952B0441C5284AA2961B1D4180BCF1FBA5E3FF401A7F4E43CDDF1A4144B16B0C1818FA40490F31E0D63718417808025F5A0FFB40690E9F2ECB04154126A94B48BB1000411D6F970D846115417EE97EED60BD0541ACF6A9CF228716413CEBA25078230C41CA8ED4B63B3C1941C50A6F700FA70E415F3769D7E8DD1A413AF85787731510411E59C611F3491D410A11BAB453351141', 1);

CREATE TABLE sub (
    cat integer,
    geom geometry(Polygon)
);

INSERT INTO sub VALUES ('010300002032BF0D0001000000070000002C456A99AD552041E490A6D20468FF406141D1B24ACC1E41AC5F98B19C06FC40C5284AA2961B1D4180BCF1FBA5E3FF409D62289E69AF1C41C89B4A27952B0441E04A333FC1311E417EE97EED60BD0541A5FC2FF222102041AE669D717D4A04412C456A99AD552041E490A6D20468FF40', 2);
INSERT INTO sub VALUES ('010300002032BF0D000100000005000000EC499BDF68082041D8A65F4924731141BB714D22A7FA1E41F7A2C0B474F31441A5FC2FF2221020411E69E2B8A15F1541CEAF8DDD62B1214176AC2241EB581441EC499BDF68082041D8A65F4924731141', 3);
INSERT INTO sub VALUES ('010300002032BF0D00010000000A000000D4012A11645C214143B25606B1C20F41886222F01CB92141D6B92362113E10417CBEE9881A6322418F515E27F8650F41FEB487FCA3FD224121C9B9B34D0D0B419E326F5946242341CC75C2F5BA950641766C4D5519B82241E1D0F7DCAC0C044114FDF8CAA8A921419531F0BB6569044128E0094DBFDF214164B4D13749DC0541E2929E5F79E72141EC7123610B160A41D4012A11645C214143B25606B1C20F41', 4);
INSERT INTO sub VALUES ('010300002032BF0D0001000000090000005F3769D7E8DD1A413AF8578773151041CA8ED4B63B3C1941C50A6F700FA70E41C6A06E07725C1A419FD21B40C4720A414F61CD13D8E2174127937A4C2AF90741ACF6A9CF228716413CEBA25078230C41789F137DE08F1541F62DA9323DFF10418A0A009A1B851841CC0227E2211D1241598DE115FFF719410CFEF59B666A12415F3769D7E8DD1A413AF8578773151041', 5);
Kompozice polygonů

Výběr překrývajících se polygonů:

 select s.cat from sub as s join ref as r on ST_Overlaps(s.geom, r.geom);

nebo

 select s.cat from sub as s join ref as r on ST_Relate(s.geom,r.geom, '2********');
 cat 
-----
   5

Procedury

Příklady uživatelských procedur pro PostGIS najdete např. zde.

Příklad funkce pro výpočet výměry v km2.

CREATE OR REPLACE FUNCTION area_km(geometry)
 RETURNS integer AS 
$$ 
SELECT ROUND(ST_Area($1)/1e6)::integer
$$
LANGUAGE 'sql';

Aplikace funkce area_km.

SELECT nazev, area_km(the_geom) FROM gis1.obce limit 3;
  nazev   | area_km 
----------+---------
 Abertamy |       9
 Adamov   |       1
 Adamov   |       3
(3 rows)

Další příklad je převzat z 3. cvičení GIS1.

CREATE OR REPLACE FUNCTION boundary_split(geometry, integer)
 RETURNS SETOF geometry AS 
$$ 
SELECT ST_Line_Substring(the_geom, $2 * n / length,
  CASE
	WHEN $2 * (n + 1) < length THEN $2 * (n + 1) / length
	ELSE 1
  END) AS the_geom
FROM 
(SELECT ST_LineMerge($1) AS the_geom,
  ST_Length($1) AS length) AS t
CROSS JOIN generate_series(0, $2) AS n
WHERE n * $2 / length < 1;
$$
LANGUAGE 'sql';

Počet liniových segmentů pro délku linie 2km.

CREATE TABLE cr AS SELECT ST_Boundary(ST_Union(the_geom)) AS the_geom FROM gis1.obce GROUP BY pomoc;

SELECT COUNT(*) FROM (SELECT boundary_split(the_geom, 2000) FROM cr) AS cr_2km;
1059

Počet liniových segmentů pro délku linie 5km.

SELECT COUNT(*) FROM (SELECT boundary_split(the_geom, 5000) FROM cr) AS cr_5km;
424

Podívejte se také na:

Související články

Externí odkazy