PostGIS
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 '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));
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"]]');
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
- Instalace PostGIS
- Cvičná databáze PostGIS
- Psql a Emacs
- PostGIS Raster
- SpatiaLite
- OpenStreetMap
- pgRouting
- PostGIS Topology
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