PostGIS: Porovnání verzí

Z GeoWikiCZ
Řá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://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://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}}
{{Databáze}}
{{GIS}}
{{GIS}}
{{GFOSS}}
{{GFOSS}}

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

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

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