Wprowadzenie do PostGIS

Przykłady geometrii WKT

• POINT(0 0)
• LINESTRING(0 0,1 1,1 2)
• POLYGON((0 0,4 0,4 4,0 4,0 0),(1 1, 2 1, 2 2,

1 2,1 1))

• MULTIPOINT(0 0,1 2)
• MULTILINESTRING((0 0,1 1,1 2),(2 3,3 2,5 4))
• MULTIPOLYGON( ((0 0,4 0,4 4,0 4,0 0),(1 1,2

1,2 2,1 2,1 1)), ((-1 -1,-1 -2,-2 -2,-2 -1,-1 -1)) )

• GEOMETRYCOLLECTION(POINT(2

3),LINESTRING((2 3,3 4)))

Informacje podstawowe

• Wymagania specyfikacji OGC:

– obiekty przestrzenne powinny być przechowywane

razem z identyfikatorem układu odniesień

przestrzennych (spatial referencing system identifier,

SRID)

– SRID jest wymagany przy tworzeniu obiektów

przestrzennych wstawianych do bazy danych

• Transformacje formatów reprezentacji obiektów

przestrzennych

– bytea WKB = asBinary(geometry);
– text WKT = asText(geometry);
– geometry = GeomFromWKB(bytea WKB, SRID);
– geometry = GeometryFromText(text WKT, SRID);

INSERT INTO SPATIALTABLE (THE_GEOM, THE_NAME )
VALUES ( GeomFromText(’POINT(-126.4 45.32)’, 312), ’A Place’ )

Informacje podstawowe

• PostGIS EWKB, EWKT …

– Specyfikacja OGC definiuje tylko geometrię 2d, a

SRID nie jest włączany w reprezentacjach we/wy

– PostGIS to rozszerza formaty OGC,

wprowadzając EWKB, EWKT (każda geometria w

postaci WKB/WKT jest poprawnym wyrażeniem

EWKB/EWKT)

– rozszerzenie może przysporzyć problemy w

przyszłości

• PostgisEWKB/EWKT

– uwzględnia współrzędne 3dm, 3dz, 4d oraz

dopuszcza włączanie informacji SRID

OGC Simple Features

Specification for SQL

• Definiuje

– standardowe typy obiektów

przestrzennych

– funkcje do manipulacji obiektami

przestrzennymi

– zbiór tabel z metadanymi

• Istnieją dwie tabele metadanych OGC

– SPATIAL_REF_SYS
– GEOMETRY_COLUMNS table

Tabela geometry_columns

Column | Type | Modifiers
--------------------+------------------------+-----------
f_table_catalog | character varying(256) | not null
f_table_schema | character varying(256) | not null
f_table_name | character varying(256) | not null
f_geometry_column | character varying(256) | not null
coord_dimension | integer | not null
srid | integer | not null
type | character varying(30) | not null
Indexes:
"geometry_columns_pk" PRIMARY KEY, btree

(f_table_catalog, f_table_schema, f_table_name,

f_geometry_column)

Tabela geometry_columns

• F_TABLE_NAME

– W pełni kwalifikowana nazwa tabeli obiektów zawierającą kolumnę z

geometrią

• F_TABLE_CATALOG, F_TABLE_SCHEMA

– „catalog” i „schema” to terminy Oracle’owe.

– „catalog” w PostgreSQL nie występuje, stąd kolumna o tej nazwie jest pusta

– „schema” w PostgreSQL występuje, stąd kolumna zawiera nazwę użytego

schematu (domyślnie public).

• F_GEOMETRY_COLUMN

– nazwa kolumny z geometrią w tablicy obiektów

• COORD_DIMENSION

– wymiar przestrzenny kolumny (2, 3 lub 4)

• SRID

– Identyfikator układu odniesień przestrzennych używanego dla współrzędnych

przechowywanej w tabeli geometrii (klucz obcy do SPATIAL_REF_SYS)

• TYPE

– typ obiektu przestrzennego przechowywanego w kolumnie z geometrią,

może być jednym z: POINT, LINESTRING, POLYGON, MULTIPOINT,

MULTILINESTRING, MULTIPOLYGON, GEOMETRYCOLLECTION

– dla kolekcji heterogenicznych można GEOMETRY jako deklaracji typu

– w nowszych wydaniach PostGIS jest jeszcze typ GEOGRAPHY

Tabela spatial_ref_sys

Column | Type | Modifiers
------------+------------------------+-----------
srid | integer | not null
auth_name | character varying(256) |
auth_srid | integer |
srtext | character varying(2048)|
proj4text | character varying(2048)|
Indexes:
"spatial_ref_sys_pkey" PRIMARY KEY, btree (srid)

SRID

• An integer value that uniquely identifies the Spatial Referencing System

(SRS) within the database. e.g., SRID = 4326 is the WGS84

AUTH_NAME

• The name of the standard or standards body that is being cited for this

reference system. For example, "EPSG" would be a valid AUTH_NAME.

AUTH_SRID

• The ID of the Spatial Reference System as defined by the Authority cited in

the AUTH_NAME. In the case of EPSG, this is where the EPSG projection

code would go.e.g., the 4326 for the “EPSG:4326” (WGS84)

SRTEXT

• The Well-Known Text representation of the Spatial Reference System.

Przykładowy SRTEXT

PROJCS["NAD83 / UTM Zone 10N",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.257222101]
],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]
],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-123],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1]
]

PROJ4TEXT

• PostGIS wykorzystuje bibliotekę Proj4

do transformacji układów odniesień

przestrzennych

• Kolumna PROJ4TEXT zawiera ciąg

znaków definiujący dany SRID

zgodnie z Proj4

– +proj=utm +zone=10 +ellps=clrk66

+datum=NAD27 +units=m

– +proj=longlat +ellps=WGS84

+datum=WGS84 +no_defs

Tworzenie tabel z danymi

przestrzennymi

1. Utwórz tabelę w „zwykły” sposób,

nazywając atrybuty wg potrzeb.

2. Dodaj kolumnę z geometrią metodą

AddGeometryColumn należącą do
rozszerzenia przestrzennego
PostGIS/OGC

3. Wstaw dane przestrzenne stosując

wyrażenia SQL insert & select

Krok 1

• Utworzenie tabeli

CREATE TABLE "public"."dublin_historical“
(gid serial PRIMARY KEY,
"rmp_prop" int8,
"map_symbol" int8,
"entity_id" varchar(7),
"co_id" int8,
"smr_val0" numeric,
"nat_grid_e" numeric,
"class_desc" varchar(255),
"nat_grid_n" numeric,
"objectid" int8,
"townlands" varchar(255),
"scope_n1" varchar(255),
"smrs" varchar(255));

W odpowiedzi PostgreSQL/PostGIS
napisze:
NOTICE: CREATE TABLE will create
implicit sequence
"dublin_historical_gid_seq" for
serial column
"dublin_historical.gid"
NOTICE: CREATE TABLE / PRIMARY
KEY will create implicit index
"dublin_historical_pkey" for table
"dublin_historical"

Krok 2

• Wstawienie kolumny z geometrią
SELECT

AddGeometryColumn('public','dub
lin_historical','the_geom','299
00','POINT',2);

W przykładzie posłużono się układem
SRID = 29900 (Irish National Grid ).
Geometrią przechowywaną w kolumnie tabeli
będzie POINT.

Krok 2

• Wstawienie ograniczeń na geometrię

(wymiar 2)

ALTER TABLE dublin_historical DROP

CONSTRAINT enforce_dims_the_geom;

ALTER TABLE dublin_historical
ADD CONSTRAINT enforce_dims_the_geom CHECK

(ndims(the_geom) = 2);

Krok 2

• Wstawienie kolejnych ograniczeń na

geometrię (POINT lub NULL)

ALTER TABLE

dublin_historical

ADD CONSTRAINT

enforce_geotype_the_geom CHECK

(geometrytype(the_geom) =

‘POINT'::text OR the_geom IS

NULL);

Krok 2

• Wstawienie kolejnych ograniczeń na

geometrię (SRID)

ALTER TABLE

dublin_historical

ADD CONSTRAINT

enforce_srid_the_geom CHECK
(srid(the_geom) = 29900);

Krok 3

• Przepisanie danych do utworzonej tabeli

insert into "dublin_historical"

("rmp_prop","map_symbol", "entity_id",

"co_id","smr_val0",

"nat_grid_e“ ,"class_desc","nat_grid_n","obj

ectid", "townlands","scope_n1", "smrs",

the_geom)

select "rmp_prop", "map_symbol", "entity_id",

"co_id","smr_val0", "nat_grid_e",

"class_desc", "nat_grid_n", "objectid",

"townlands", "scope_n1", "smrs", h.the_geom

from "county" c, historical h where

contains(c.the_geom,h.the_geom) and

name = 'Dublin';

Załadowanie danych z pliku i

eksport do pliku

• użycie shp2pgsql (konwersja ESRI Shape do SQL)

# shp2pgsql shaperoads myschema.roadstable >

roads.sql

• użycie komendy psql:

# psql -d roadsdb -f roads.sql

• użycie shp2pgsql (łączącego się z bazą danych,

konwertującego tabelę, zdefiniowaną opcjonalnym

zapytaniem, do pliku shp)

# pgsql2shp [<options>] <database>

[<schema>.]<table>

# pgsql2shp [<options>] <database> <query>
(<options>: host, user, password, filename, port,

geometry column, etc.)

Operatory

&& : sprawdza, czy BBOX jednej

geometrii przecina się z BBOX drugiej

~= : sprawdza, czy dwie geometrie są

identyczne geometrycznie

= : sprawdzam czy BBOXy dwóch

geometrii są takie same

Indeksy wspierane przez

PostGIS

• B-Trees

– They are used for data which can be sorted along one axis (no

use for GIS data)

• R-Trees

– break up data into rectangles, and sub-rectangles, and sub-sub

rectangles, etc. R-Trees are used by some spatial databases to

index GIS data, but the PostgreSQL R-Tree implementation is

not as robust as the GiST implementation

• GiST (Generalized Search Trees)

– break up data into "things to one side", "things which overlap",

"things which are inside" and can be used on a wide range of

data-types, including GIS data. PostGIS uses an R-Tree index

implemented on top of GiST to index GIS data

CREATE INDEX [indexname] ON [tablename]
USING GIST ( [geometryfield] GIST_GEOMETRY_OPS );

GiST (Generalized Search

Trees)

• only the bounding-box-based operators such as

&& can take advantage of the GiST spatial index

– This is very slow…
SELECT the_geom FROM geom_table
WHERE distance( the_geom,

GeomFromText( ’POINT(100000

200000)’, -1 ) ) < 100

– If a GiST is available over the_geom, the following is

much more efficient...

SELECT the_geom FROM geom_table
WHERE the_geom && ’BOX3D(90900 190900, 100100
200100)’::box3d
AND distance( the_geom,

GeomFromText( ’POINT(100000

200000)’, -1 ) ) < 100

PostGIS (OGC) Functions

• Management Functions

– AddGeometryColumn(varchar, varchar,

varchar, integer,varchar, integer)

– DropGeometryColumn(varchar, varchar,

varchar))

– SetSRID(geometry)

PostGIS (OGC) Functions

• Geometry Relationship Functions

– Distance(geometry,geometry)

– Equals(geometry,geometry)

– Disjoint(geometry,geometry)

– Intersects(geometry,geometry)

– Touches(geometry,geometry)

– Crosses(geometry,geometry)

– Within(geometry,geometry)

– Overlaps(geometry,geometry)

– Contains(geometry,geometry)

– Relate(geometry,geometry,

intersectionPatternMatrix)

– Relate(geometry,geometry)

PostGIS (OGC) Functions

• Geometry Processing Functions

– Centroid(geometry)

– Area(geometry)

– Length(geometry)

– PointOnSurface(geometry)

– Boundary(geometry)

– Buffer(geometry,double,[integer])

– ConvexHull(geometry)

– Intersection(geometry,geometry)

– SymDifference(geometry,geometry)

– Difference(geometry,geometry)

– GeomUnion(geometry,geometry)

– GeomUnion(geometry set)

– MemGeomUnion(geometry set)

PostGIS (OGC) Functions

• Geometry Accessors

– AsText(geometry)

– AsBinary(geometry)

– SRID(geometry)

– Dimension(geometry)

– Envelope(geometry)

– IsEmpty(geometry)

– IsSimple(geometry)

– IsClosed(geometry)

– IsRing(geometry)

– NumGeometries(geometry)

– GeometryN(geometry,int)

PostGIS (OGC) Functions

• Geometry Accessors

– NumPoints(geometry)

– PointN(geometry,integer)

– ExteriorRing(geometry)

– NumInteriorRings(geometry)

– InteriorRingN(geometry,integer)

– EndPoint(geometry)

– StartPoint(geometry)

– GeometryType(geometry)

– X(geometry)

– Y(geometry)

– Z(geometry)

PostGIS (OGC) Functions

• Geometry Constructors

– GeomFromText(text,[<srid>])
– PointFromText(text,[<srid>])
– LineFromText(text,[<srid>])
– LinestringFromText(text,[<srid>])
– PolyFromText(text,[<srid>])
– PolygonFromText(text,[<srid>])
– MPointFromText(text,[<srid>])
– MLineFromText(text,[<srid>])
– MPolyFromText(text,[<srid>])
– GeomCollFromText(text,[<srid>])

PostGIS (OGC) Functions

• Geometry Constructors (cont’)

– GeomFromWKB(bytea,[<srid>])
– PointFromWKB(bytea,[<srid>])
– LineFromWKB(bytea,[<srid>])
– LinestringFromWKB(bytea,[<srid>])
– PolyFromWKB(bytea,[<srid>])
– PolygonFromWKB(bytea,[<srid>])
– MPointFromWKB(bytea,[<srid>])
– MLineFromWKB(bytea,[<srid>])
– MPolyFromWKB(bytea,[<srid>])
– GeomCollFromWKB(bytea,[<srid>])

PostGIS Extensions

• PostGIS also offers … functions not

supported by OGC

– Management Functions
– Operators
– Measurement Functions
– Geometry Outputs
– Geometry Constructors
– Geometry Editors
– Misc

Sprawdzanie poprawności

geometrii

select IsValid(the_geom) from lakes;
NOTICE: Self-intersection at or near point 114275 271699
NOTICE: Self-intersection at or near point 124552 240642
NOTICE: Self-intersection at or near point 121283 305664
isvalid
---------
t
t
f
f
f
t
t
t

Punkt centralny hrabstwa

select name,

asText(Centroid(the_geom)) from
county;

Znajdowanie prostokąta

ograniczającego

• BBOX dla każdego hrabstwa:

select name,

asText(envelope(the_geom)) from county;

lub

select name, extent(the_geom) from

county group by name;

• group by pozwala na przetworzenie wielu

poligonów, extent zwraca geometrię jako

tekst.

Znajdowanie wymiaru

prostokąta ograniczającego

• Na początek można znaleźć BBOX

select name, extent(the_geom) from county

where name= 'Dublin' group by name;

• wynik:

Dublin | BOX(297049.5 215877.484375,330201.8125

266854.40625)

• mając współrzędne narożników można obliczyć ich

odległość:

select distance(
geomFromText(‘Point….)’),
geomFromText(‘Point….)’));

Znajdowanie punktów w

wielokącie

• Do tego zadania można użyć

niestandardowego operatora &&.

• SELECT name FROM original_county

where GeomFromText('POINT(309612.0

233192.0)', 29900) && the_geom;

• Tworząc zapytania należy pamiętać, że

tylko operatory bazujące na BBOX jak && w

pełni mogą wykorzystać indeks

przestrzenny GiST PostgreSQL’a

Znajdowanie obiektów w

okolicach punktu

SELECT townlands FROM dublin_historical

WHERE distance(the_geom,

GeomFromText('POINT(309612.0 233192.0)',

29900)) < 100;

SELECT c.name, h.townlands
FROM county AS c, dublin_historical AS h

WHERE distance(h.the_geom,

GeomFromText('POINT(317431.0 231704.0)',

29900)) < 1000 and c. Name = 'Dublin County

Borough';

Znajdowanie obiektów w

okolicach punktu

SELECT c.name, h.townlands
FROM county AS c,dublin_historical AS h
WHERE st_dwithin(h.the_geom,

PointFromText('POINT(317431.0 231704.0)',

29900),1000) and c. Name = '

Dublin County

Borough

';

SELECT gid,townlands
FROM dublin_historical
WHERE st_dwithin(the_geom,

GeomFromText('POINT(317431.0 231704.0)',29900),500)

ORDER BY st_dwithin(the_geom,

GeomFromText('POINT(317431.0 231704.0)',29900),500)

LIMIT 100;

Znajdowanie największego

hrabstwa

SELECT name, area(the_geom)/10000 AS

hectares FROM county ORDER BY hectares

DESC LIMIT 1;

SELECT name, area(the_geom)/10000 AS

hectares FROM county

WHERE name != 'Northern Ireland‘ ORDER BY

hectares DESC LIMIT 1;

Analizy przestrzenne z drogami

• Długość dróg w pełni zawartych w Dublin Borough
SELECT c.name, sum(length(r.the_geom))/1000

as roads_km FROM roads AS r, county AS c

WHERE r.the_geom && r.the_geom

AND contains(c.the_geom,r.the_geom)
AND c.name = 'Dublin County Borough'
GROUP BY c.name
ORDER BY roads_km;

• Obiekty historyczne w pobliżu dróg hrabstwa Dublin
SELECT townlands FROM dublin_historical h,

roads r, county c WHERE

distance(h.the_geom,r.the_geom) < 1000
AND c.name = 'Dublin County Borough';

Nakładanie się obiektów

• Możliwe są testy sprawdzające

przecinanie się obiektów oraz
operacje przecięcia:

ST_Intersects(a,b)=TRUE | FALSE
ST_Intersection(a,b)=GEOMETRY

Niuanse analiz

przestrzennych

Niuanse analiz przestrzennych

• Oryginalne zapytanie:
SELECT the_geom, name
FROM neighborhoods
WHERE name IN('Hyde Park','Roxbury')

• expand: tworzy prostokąt, który rozszerza o 500 jednostek bbox dla każdej

geometrii

SELECT ST_Expand(the_geom,500) as geom, name
FROM neighborhoods
WHERE name IN('Hyde Park','Roxbury')

• extent: tworzy pojedynczy bbox wokół wybranych geometrii
SELECT ST_Extent(the_geom) as thebbox, name
FROM neighborhoods
WHERE name IN('Hyde Park','Roxbury')

• buffer: rozszerza każdą geometrię o 500 jednostek
SELECT ST_Buffer(the_geom,500) as geom, name
FROM neighborhoods
WHERE name IN('Hyde Park','Roxbury')

Expand

--box meter expanded box around bbox of a linestring
SELECT CAST(ST_Expand(ST_GeomFromText(’LINESTRING(2312980

110676,2312923 110701,2312892 110714)’, 2163),10) As box2d);

st_expand
------------------------------------
BOX(2312882 110666,2312990 110724)
--10 meter expanded 3d box of a 3d box
SELECT ST_Expand(CAST(’BOX3D(778783 2951741 1,794875

2970042.61545891 10)’ As box3d),10)

st_expand
-----------------------------------------------------
BOX3D(778773 2951731 -9,794885 2970052.61545891 20)
--10 meter geometry astext rep of a expand box around a point geometry
SELECT

ST_AsEWKT(ST_Expand(ST_GeomFromEWKT(’SRID=2163;POINT(2312980

110676)’),10));

st_asewkt
------------------------------------------------------------------------------------------------- -
SRID=2163;POLYGON((2312970 110666,2312970 110686,2312990

110686,2312990 110666,2312970 -

110666))

Extent

SELECT ST_Extent(the_geom) as bextent FROM sometable;
st_bextent
------------------------------------
BOX(739651.875 2908247.25,794875.8125 2970042.75)
--Return extent of each category of geometries
SELECT ST_Extent(the_geom) as bextent
FROM sometable
GROUP BY category ORDER BY category;
bextent | name
----------------------------------------------------+----------------
BOX(778783.5625 2951741.25,794875.8125 2970042.75) | A
BOX(751315.8125 2919164.75,765202.6875 2935417.25) | B
BOX(739651.875 2917394.75,756688.375 2935866) | C
--Force back into a geometry
PostGIS 1.5.2 Manual
262 / 317
-- and render the extended text representation of that geometry
SELECT ST_SetSRID(ST_Extent(the_geom),2249) as bextent FROM sometable;
bextent
--------------------------------------------------------------------------------
SRID=2249;POLYGON((739651.875 2908247.25,739651.875 2970042.75,794875.8125

2970042.75,

794875.8125 2908247.25,739651.875 2908247.25))

Buffer

--A buffered point approximates a circle
SELECT ST_NPoints(ST_Buffer(ST_GeomFromText(’POINT(100 90)’), 50))

As promisingcircle_pcount,

ST_NPoints(ST_Buffer(ST_GeomFromText(’POINT(100 90)’), 50, 2)) As

lamecircle_pcount;

promisingcircle_pcount | lamecircle_pcount
------------------------+-------------------
33 | 9

Zapytania o sieć drogową

Najkrótsza droga z

Dublina do

Waterford

• Możliwe dzięki pgRouting

Przegląd dostępnych w

pgRouting zapytań za pomocą

pgAdmin III

Network Queries

-- The start node 381 is in Dublin
-- The destination is node 660 in Waterford
-- Must return columns 'id', 'source', 'target' and
'cost'
SELECT * FROM shortest_path('
SELECT gid as id,
source::integer,
target::integer,
length::double precision as cost
FROM all_roads_e',
381, 660, false, false);

Network Queries

The edges are stored in all_roads_e
The source and target information was created in
OpenJump using Tools | Analysis | Planar Graph.
Details will be described in labs:

----------+------------------+--------------------
gid | integer | not null default nextval('egclass)
id | bigint
source | bigint
target | bigint
the_geom | geometry
length | double precision

Długość wszystkich odcinków

dróg w pełni zawartych w

hrabstwie

SELECT  

  Sum(ST_Length(r.the_geom))/1000 

    AS kilometers  

FROM  

  all_roads r,  

dublin_regions dr 

WHERE  

  ST_Contains(dr.the_geom,

r.the_geom) AND 

   dr.name = ‘Dublin County

Borough';

Zapytania

• Dwa domy w odległości 500 metrów od biblioteki z

największą liczbą mieszkańców

SELECT b1.residents
FROM buildings b1, buildings b2
WHERE
b2.name = 'CHESTER BEATTY LIBRARY' and
ST_DWithin(b2.the_geom,b1.the_geom, 500) ORDER BY

residents DESC LIMIT 2;

• Hrabstwa mające dokładnie jednego sąsiada
SELECT c1.name
FROM all_counties c1, all_counties c2
WHERE touches(c1.the_geom, c2.the_geom) = 'TRUE'
GROUP BY c1.name
HAVING count(c2.name) =1;

Szerokość geograficzna

budynków wysuniętych

najbardziej na północ i na

południe

SELECT st_y(st_transform(the_geom,4326)) AS

latitude from buildings ORDER BY latitude

DESC limit 1;

SELECT st_y(st_transform(the_geom,4326)) AS

latitude from buildings ORDER BY latitude

ASC limit 1;

SELECT st_y(the_geom) AS easting FROM

buildings ORDER BY easting DESC limit 1;

Zapytania arytmetyczne

SELECT saps_label,primary_degree10_4,(male1_1 + female1_1) as Pop,
((cast (primary_degree10_4 as float)) / (cast ((male1_1 +

female1_1) as float)) * 100) as Grad_percent

FROM dublin_eds
WHERE primary_degree10_4 IS NOT NULL
ORDER BY Grad_percent
DESC LIMIT 1;

saps_label | primary_degree10_4 | pop | grad_percent
--------------------+--------------------+------+--------------
117 Mansion House A | 279 | 3802 |7.33824302998422