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 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#Poznámky.
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
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