PostGIS: Porovnání verzí

Z GeoWikiCZ
m (Obsah stránky nahrazen textem „{{freegiswiki|PostGIS}}“)
 
Řádek 1: Řádek 1:
[[Image:postgis.png|150px|right]]
{{freegiswiki|PostGIS}}
'''[http://www.postgis.org PostGIS]''' je rozšíření objektově-relačního databázového systému [[PostgreSQL]] pro podporu geografických objektů. PostGIS implementuje specifikaci [http://www.opengeospatial.org/standards/sfa Simple Features] konsorcia [http://www.opengeospatial.org Open Geospatial Consortium].
 
;Rozšíření:
 
* [[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čná databáze PostGIS|cvičnou databázi]]''' na serveru 'geo102'.
 
== Příklady ==
 
=== [[153GIS1|GIS 1]] ===
 
* [http://geo.fsv.cvut.cz/~gin/uzpd/examples/postgis/gis1-2.sql 2. cvičení]
* [http://geo.fsv.cvut.cz/~gin/uzpd/examples/postgis/gis1-3.sql 3. cvičení]
* [http://geo.fsv.cvut.cz/~gin/uzpd/examples/postgis/gis1-4.sql 4. cvičení]
 
=== [[153UZPD|Ú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:
 
<source lang=bash>
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
</source>
 
'''Poznámka:''' Pro nahrání dávky <tt>postgis.sql</tt> jsou vyžadována práva "superuživatele" (superuser) databázového systému [[PostgreSQL#Poznámky|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í <code>PGCLIENTENCODING</code>, např.
 
<source lang="bash">
export PGCLIENTENCODING=WIN1250
</source>
 
=== SQL ===
 
<source lang="sql">
INSERT INTO zb (cat, zb_geom) VALUES (1, ST_GeomFromText('POINT(13.30150173 49.79076912)', 4326));
</source>
 
nebo
 
<source lang="sql">
INSERT INTO zb (cat, zb_geom) VALUES (1, ST_SetSRID(ST_MakePoint(13.30150173,49.79076912), 4326));
</source>
 
Viz [http://www.postgis.org/documentation/manual-svn/ch07.html#Geometry_Constructors 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 &rarr; 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) ===
 
{{GrassPrikaz|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'.
 
<source lang='python'>
#!/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())
</source>
 
=== Další formáty podporované knihovnou [[GDAL/OGR|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 <tt>geometry_columns</tt>, 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ř.
<source lang="sql">
SELECT AddGeometryColumn('zb', 'zb_geom', 4326, 'POINT', 2);
</source>
 
*;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ř.
 
<source lang="sql">
SELECT DropGeometryTable('osm', 'czech_roads');
</source>
 
Při manuálním vytvoření prostorového atributu (např. <code>CREATE TABLE tabulka AS ...</code>) je potřeba aktualizovat tabulku <tt>geometry_columns</tt> (INSERT nebo pomocí funkce <tt>Populate_Geometry_Columns</tt> (od verze 1.4.0)). Např.:
 
<source lang="sql">
SELECT Populate_Geometry_Columns('public.zb'::regclass);
</source>
 
=== Definice S-JTSK (ESRI:102067) ===
 
<source lang="sql">
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.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56',
'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]]');
</source>
 
Viz [http://spatialreference.org/ref/esri/102067/ SpatialReference.org]
 
=== Rozsah dat ===
 
<source lang="sql">
SELECT ST_Extent(the_geom) from points;
-- EPSG:4326 --
SELECT ST_Extent(ST_Transform(the_geom, 4326)) from points;
</source>
 
=== Self-intersection ===
 
Příklad detekce průniku polygonů lesních porostů (data [[OpenStreetMap]]).
 
<source lang="sql">
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));
</source>
 
=== Round() ===
 
Jaká je výměra lesa na mapovém listu 'M-33-63-A' v km2?
 
<source lang=sql>
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);
</source>
 
125
 
<source lang=sql>
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);
</source>
 
124.71
 
=== Prvek &rarr; multi-prvek ===
 
Příklad pro převod 'point' na 'multipoint' na základě kategorie (<tt>cat</tt>). Vstupní vrstva
 
<pre>
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)
</pre>
 
<source lang="sql">
SELECT cat,st_astext(st_union(wkb_geometry)) from points group by cat;
</source>
 
<pre>
cat | st_numgeometries
-----+------------------
  1 |                7
  2 |                3
(2 rows)
</pre>
 
----
 
Příklad pro tabulku obcí ze [[Cvičná databáze PostGIS|cvičné databáze PostGIS]].
 
<source lang=sql>
SET search_path TO public,gis1;
</source>
 
Počet obcí v ČR
 
<source lang=sql>
SELECT count(DISTINCT kodob) from obce;
</source>
<pre>
count
-------
  6249
</pre>
 
Řada obcí je tvořena více polygony
 
<source lang=sql>
SELECT nazev,count(*) from obce GROUP BY kodob,nazev HAVING count(*) > 1;
</source>
<pre>
        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)
</pre>
 
Cílem je vytvořit novou tabulku obcí, kde pro tyto záznamy bude místo 'polygonu' použit 'multipolygon', tj. bude platit "<tt>jeden záznam == jedna obec</tt>".
 
<source lang=sql>
-- 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);
</source>
 
=== 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ř. [http://trac.osgeo.org/postgis/wiki/UsersWikiplpgsqlfunctions zde].
 
Příklad funkce pro výpočet výměry v km<sup>2</sup>.
 
<source lang="sql">
CREATE OR REPLACE FUNCTION area_km(geometry)
RETURNS integer AS
$$
SELECT ROUND(ST_Area($1)/1e6)::integer
$$
LANGUAGE 'sql';
</source>
 
Aplikace funkce <tt>area_km</tt>.
 
<source lang="sql">
SELECT nazev, area_km(the_geom) FROM gis1.obce limit 3;
</source>
 
<pre>
  nazev  | area_km
----------+---------
Abertamy |      9
Adamov  |      1
Adamov  |      3
(3 rows)
</pre>
 
Další příklad je převzat z [http://geo.fsv.cvut.cz/~gin/uzpd/examples/postgis/gis1-3.sql 3. cvičení GIS1].
 
<source lang="sql">
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';
</source>
 
Počet liniových segmentů pro délku linie 2km.
 
<source lang="sql">
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;
</source>
 
<pre>
1059
</pre>
 
Počet liniových segmentů pro délku linie 5km.
 
<source lang="sql">
SELECT COUNT(*) FROM (SELECT boundary_split(the_geom, 5000) FROM cr) AS cr_5km;
</source>
 
<pre>
424
</pre>
 
Podívejte se také na:
 
* [http://www.postgres.cz/index.php/PL/pgSQL PL/pgSQL]
* [http://www.root.cz/clanky/jemny-uvod-do-plpgsql/ Jemný úvod do PL/pgSQL]
* [http://www.postgresql.org/docs/8.3/static/plpgsql.html PL/pgSQL - SQL Procedural Language]
 
== Související články ==
 
* [[Instalace PostGIS]]
* [[Cvičná databáze PostGIS]]
* [[Psql a Emacs]]
* [[SpatiaLite]]
* [[OpenStreetMap]]
 
== Výuka ==
 
* [[153UZPD|Úvod do zpracování prostorových dat]]
* [[153YFSG|Free software GIS]]
 
== Externí odkazy ==
 
* [http://www.postgis.org/documentation/ PostGIS documentation]
* [http://www.foss4g2007.org/workshops/W-04/ Introduction to PostGIS] by Paul Ramsey
* [http://revenant.ca/www/postgis/workshop/index.html Introduction to PostGIS (OSGeo)] by Mark Leslie
* [http://www.mapbender.org/presentations/Spatial_Data_Management_Arnulf_Christl/Spatial_Data_Management_Arnulf_Christl.pdf Spatial Data Management]
* [http://edndoc.esri.com/arcsde/9.0/general_topics/understand_spatial_relations.htm Understanding spatial relations] by ESRI
* [http://www.bostongis.com/postgis_quickguide.bqg PostGIS Quick Guide - Cheatsheet]
* [http://grass.osgeo.org/wiki/Spatial_SQL Spatial SQL on GRASSWiki]
* [http://www.manning.com/obe/ PostGIS in action]
* [http://pgrouting.postlbs.org pgRouting]
* [http://spatialreference.org Spatial Reference]
* [http://trac.osgeo.org/postgis/wiki/UsersWikiPostgisTopology PostGIS Topology]
* [http://s3.opengeo.org/postgis-power.pdf Tips for the PostGIS Power User] (FOSS4G 2010 - Paul Ramsey)
* [http://workshops.opengeo.org/postgis-spatialdbtips/ Spatial database tips and tricks]
* [http://www.kappasys.ch/cms/index.php?id=23&L=5 PostGIS Versioning - pgVersion]
* [http://planet.postgis.net/ PostGIS Planet] (blogy)
* [http://workshops.opengeo.org/postgis-intro/index.html Introduction to PostGIS] (OpenGeo)
 
 
<!-- * [http://gis.zcu.cz/studium/ugi/referaty/05/PostGIS/index.html Představení projektu PostGIS] -->
* [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 ČVUT v Praze)
 
{{Databáze}}
{{GIS}}
{{GFOSS}}

Aktuální verze z 8. 1. 2013, 14:32

Stránky přesunuty na Free GIS Portál: http://geo.fsv.cvut.cz/freegis/PostGIS