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.

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

Příklady

153GIS1

153UZPD

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

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

Citováno z „https://geo.fsv.cvut.cz/index.php?title=PostGIS&oldid=18505“