PostGIS: Porovnání verzí

Verze z 10. 4. 2012, 11:09

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.

Nadstavby

PostGIS Raster - uložení, manipulace s rastrovými daty v prostředí PostGIS
PostGIS Topology - topologická správa vektorových dat
pgRouting - síťové analýzy v PostGISu

Pro experimenty používejte cvičnou databázi na serveru 'geo102'.

Příklady

GIS 1

Úvod do zpracování prostorových dat

http://geo.fsv.cvut.cz/~gin/uzpd/examples/postgis/

Vytvoření databáze

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

Příklad založení PostGIS databáze:

export DB=template_postgis
export PGPATH=/usr/share/postgresql/8.4/contrib/postgis-2.0/

createdb $DB
createlang plpgsql $DB
psql -d $DB -f $PGPATH/postgis.sql
psql -d $DB -f $PGPATH/spatial_ref_sys.sql 

# volitene PostGIS Raster
psql -d $DB -f $PGPATH/rtpostgis.sql

# volitene PostGIS Topology
psql -d $DB -f $PGPATH/topology.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 (v.out.ogr)

v.out.ogr input=cr type=area format=PostgreSQL dsn="pg:dbname=pgis_yfsg"

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

Atribut geometrie

Pro manipulaci s atributem geometrie slouží specializované funkce, které aktualizují tabulku geometry_columns, vytvoří nad atributem prostorový index.

Vytvoření prostorového atributu

SELECT AddGeometryColumn(
schéma,                        # volitelný atribut
název tabulky,
sloupec,
srid,
typ geometrie,
dimenze);

Např.

SELECT AddGeometryColumn('zb', 'zb_geom', 4326, 'POINT', 2);

Odstranění atributu geometrie

DropGeometryColumn(
schéma,                        # volitelný atribut
název tabulky,
sloupec)

Odstranění tabulky

DropGeometryTable(
schéma,                        # volitelný atribut
název tabulky)

Např.

SELECT DropGeometryTable('osm', 'czech_roads');

Při manuálním vytvoření prostorového atributu (např. CREATE TABLE tabulka AS ...) je potřeba aktualizovat tabulku geometry_columns (INSERT nebo pomocí funkce Populate_Geometry_Columns (od verze 1.4.0)). Např.:

SELECT Populate_Geometry_Columns('public.zb'::regclass);

Definice S-JTSK (ESRI:102067)

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values
( 102067, 'esri', 102067, '+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',
'PROJCS["krovak",GEOGCS["bessel",DATUM["unknown",SPHEROID
["Bessel_1841",6377397.155,299.1528128],TOWGS84[570.8,85.7,462.8,4.998,1.587,5.261,3.56]],
PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Krovak"],
PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",0],PARAMETER["azimuth",0],
PARAMETER["pseudo_standard_parallel_1",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],
PARAMETER["false_northing",0],UNIT["Meter",1]]');

Viz SpatialReference.org

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);

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);

Kopírovat data mezi databázemi

Příklad kopírování dat definovaných schématem mezi různými databázemi.

pg_dump -Fc -b -v -f gis1.dump -n gis1 pgis_student
pg_restore -Fc -d pgis_uzpd gis1.dump

Procedury

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

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

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;

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

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

Podívejte se také na:

Související články

Výuka

Externí odkazy

PostGIS documentation
Introduction to PostGIS by Paul Ramsey
Introduction to PostGIS (OSGeo) by Mark Leslie
Spatial Data Management
Understanding spatial relations by ESRI
PostGIS Quick Guide - Cheatsheet
Spatial SQL on GRASSWiki
PostGIS in action
pgRouting
Spatial Reference
PostGIS Topology
Tips for the PostGIS Power User (FOSS4G 2010 - Paul Ramsey)
Spatial database tips and tricks
PostGIS Versioning - pgVersion

PostGIS pro vývojáře
Správa časoprostorových dat v prostředí PostgreSQL/PostGIS
SQL-MM-SPATIAL (FEL ČVUT v Praze)

@@ Řádek 488: / Řádek 488: @@
 * [http://geoinformatics.fsv.cvut.cz/gwiki/PostGIS_pro_v%C3%BDvoj%C3%A1%C5%99e PostGIS pro vývojáře]
 * [http://www.openweekend.cz/slides/ow_2005/orlik_postgis.pdf Správa časoprostorových dat v prostředí PostgreSQL/PostGIS]
-* [http://service.felk.cvut.cz/courses/X36SQL/Prednasky/SQL-MM-SPATIAL.pdf SQL-MM-SPATIAL] (FEL CVUT v Praze)
+* [http://service.felk.cvut.cz/courses/X36SQL/Prednasky/SQL-MM-SPATIAL.pdf SQL-MM-SPATIAL] (FEL ČVUT v Praze)
 {{Databáze}}
 {{GIS}}
 {{GFOSS}}