PostGIS: Porovnání verzí
m (→Procedury) |
m (→Procedury) |
||
| Řádek 193: | Řádek 193: | ||
<source lang="sql"> | <source lang="sql"> | ||
CREATE OR REPLACE FUNCTION area_km( | CREATE OR REPLACE FUNCTION area_km(geometry) | ||
RETURNS integer AS | RETURNS integer AS | ||
$$ | $$ | ||
Verze z 7. 5. 2010, 08:55
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/lwpostgis.sql psql -d <databáze> -f cesta/k/spatial_ref_sys.sql
Příklad skriptu pro založení PostGIS databáze:
#!/bin/sh
if test -z "$1" ; then
echo "Usage: $0 database"
exit 1
fi
DB=$1
PGPATH=/usr/share/postgresql/8.1/contrib/
createdb $DB
createlang plpgsql $DB
psql -d $DB -f $PGPATH/lwpostgis.sql
psql -d $DB -f $PGPATH/spatial_ref_sys.sql
exit 0
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"
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
( 9102067, '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)
Příklad převzat z 3. cvičení GIS1.
CREATE OR REPLACE FUNCTION boundary_split(geom geometry, segment integer)
RETURNS 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';
Podívejte se také na: