PostGIS

Z GeoWikiCZ

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

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

Příklady

GIS 1

Úvod do zpracování prostorových dat

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=pgisďb
export PGPATH=/usr/share/postgresql/8.4/contrib/postgis-1.5/

createdb $DB
createlang plpgsql $DB
psql -d $DB -f $PGPATH/postgis.sql
psql -d $DB -f $PGPATH/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 (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 +lat_0=49.5 +lon_0=24.83333333333333
+alpha=30.28813975277778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs ',
'PROJCS["S-JTSK_Krovak_East_North",GEOGCS["GCS_S_JTSK",
DATUM["Jednotne_Trigonometricke_Site_Katastralni",SPHEROID["Bessel_1841",6377397.155,299.1528128]],
PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Krovak"],
PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Pseudo_Standard_Parallel_1",78.5],
PARAMETER["Scale_Factor",0.9999],PARAMETER["Azimuth",30.28813975277778],
PARAMETER["Longitude_Of_Center",24.83333333333333],PARAMETER["Latitude_Of_Center",49.5],
PARAMETER["X_Scale",-1],PARAMETER["Y_Scale",1],PARAMETER["XY_Plane_Rotation",90],UNIT["Meter",1],
AUTHORITY["EPSG","102067"]]');

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

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

Výuka

Externí odkazy