- 81 Views
- Uploaded on
- Presentation posted in: General

EnterpriseDB Corporation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

EnterpriseDB Corporation

EnterpriseDB is the leading provider of enterprise-class products and services based on PostgreSQL, the world's most advanced open source database.

Over 600 customers including: Sony, FTD, British Telecom, TDAmeritrade

Introduction to PostGIS

Installation – Tutorial - Exercises

What is PostGIS?

- Like Oracle Spatial, DB2 Spatial, and SQL Server Spatial, PostGIS adds spatial capabilities to an existing relational database, in this case, PostgreSQL.
- It could be renamed to “PostgreSQL Spatial” as it functions in the same way as proprietary database extensions.

What is PostGIS?

- Adds support for geographic objects to the PostgreSQL object-relational database
- PostgreSQL already has “geometric types” but native geometries are too limited for GIS data and analysis

What is PostGIS?

- PostGIS adds an indexing mechanism to allow queries with spatial restrictions or filters (“within this bounding box”) to return records very quickly from large tables.

What is PostGIS?

- PostGIS adds a “geometry” data type to the usual data types (ie. varchar, integer, date, etc.)
- PostGIS adds spatial predicates and functions that use the geometry data type.
- ST_Distance(geometry, geometry)
- ST_Area(geometry)
- ST_Intersects(geometry, geometry))

PostGIS Overview / History

- Open source
- General Public License (GPL)
- Open development and support

- History
- 2001: First release, Mapserver support
- 2002: Improved functions, indexes
- 2003: GEOS support, many functions
- 2004: SFSQL conformance
- 2005: Lightweight geometries
- 2006: OpenGIS SFSQL compliance
- 2007: SQL/MM, curves & performance
- 2008-2010: Performance enhancements

Why PostGIS?

- Because databases are better than files
- Unified Storage, Management, Access
- SQL Everywhere

- Transactional Integrity
- Multiple Users, Multiple Edits

PostGIS

PostGIS in the Spatial Stack

LAN

Internet

Mapserver

uDig

OpenIMF

GeoServer

GRASS

WebClient

QGIS

MapGuide

ArcGIS

OpenJUMP

uDig

Who is using PostGIS?

- Lots of people …
- 1790 mailing list members
- 14K visits/month10K visitors/month
- 130 source code downloads per day
- 970K Google search results
- Great Google trends …

Who is using PostGIS?

- North Dakota State Water Commission in Bismark uses PostGIS as the cost-effective foundation for managing the use of state-wide water resources. PostGIS tied unobtrusively into their existing desktop systems and pre-empted a costly migration to an ArcIMS/ArcSDE system.

Who is using PostGIS?

- The Ministry of Sustainable Resource Management (British Columbia, Canada) uses PostGIS to store, manage and analyze their Digital Road Atlas, a large and complex road network database.

Who is using PostGIS?

- Institut Géographique National (IGN France) uses PostGIS to underpin a system containing more than 100 million spatial objects, including a 3D national database of roads, rails, hydrography, vegetation, buildings and administrative boundaries.

Who is using PostGIS?

- InfoTerra (United Kingdom) uses PostGIS to store >600 million features of topographic data. Sits at the back-end of their online DataStore. Notable load times of 12 hours for the entire data-suite ~14000 features per second. More notable savings of ~£1000000 per annum.

Workshop Summary

- Introduction to Spatial Concepts
- Installation (PostgreSQL and PostGIS)
- Creating and Loading Data
- Creating Indexes
- Spatial Operations and Predicates
- Spatial SQL
- System Configuration

- A Geographical Information System (GIS) is any system used for capturing, storing, analyzing, managing, or presenting geospatial data.

- Geometric Datatypes
- Geometry Validity

2.1 Geometric Datatypes

- Each geometry has a Well-Known Text representation (WKT)
- A geometry type
- A comma-separated list of coordinate pairs

2.1 Geometric Datatypes

- Some examples are:
POINT ( 7 6 )

2.1 Geometric Datatypes

- Some examples are:
MULTIPOINT ( 1 4, 1 5, -0.5 5.5314 )

2.1 Geometric Datatypes

- Some examples are:
LINESTRING ( 1 3, 1 1, 3 0 )

2.1 Geometric Datatypes

- Some examples are:
MULTILINESTRING (( 2 1, 2 2, 3 4 ),

( 4 3, 3 2, 4 1, 3 3 ))

2.1 Geometric Datatypes

- Some examples are:
POLYGON (( 2 5, 3 8, 6 8, 5 4, 2 5 ),

( 3 6, 4 5, 5 7, 4 6, 3 6 ))

2.1 Geometric Datatypes

- Some examples are:
MULTIPOLYGON ((( 6 2, 6 4, 8 4, 9 1, 6 2 ),

( 7 3, 8 2, 8 3, 7 3 )),

(( 9 2, 8 5, 10 4, 10 3, 9 2 )))

2.1 - Exercises

- Time for some exercises!

2.3 - Geometry Validity

- PostGIS is compliant with the Open Geospatial Consortium’s (OGC) OpenGIS Specifications
- Most functions require / assume geometries are valid and in many cases, simple.
- Geometry validity really only applies for areal geometries
- Geometry simplicity applies to point and linear geometies

2.3 - Geometry Validity

- A POINT is inherently simple as a 0-dimentional geometry object.
- A MULTIPOINT is simple if no two POINTs are the same.

2.3 - Geometry Validity

- By definition, a LINESTRING is always valid
- It is simple if it does not pass through the same point twice.

SIMPLE

NOT

SIMPLE

SIMPLE

NOT

SIMPLE

2.3 - Geometry Validity

- MULTILINESTRINGs are simple if all its elements are simple and they only intersect at boundary points.

NOT SIMPLE

SIMPLE

SIMPLE

2.3 - Geometry Validity

- By definition, a POLYGON is always simple.
- It is valid if
- The boundary is made up of simpleLINEARRINGs
- boundary rings don’t cross
- holes are completely within the outer ring and touch the outer ring at most one point.
- it does not contain cutlines / spikes

2.3 - Geometry Validity

- These POLYGONs are …

valid

invalid

invalid

invalid

valid

invalid

2.3 - Geometry Validity

- By definition, a MULTIPOLYGON is valid iff
- All elements are valid
- Element interiors cannot intersect
- Element boundaries can touch, but only at a finite number of POINTs.

2.3 - Geometry Validity

- These MULTIPOLYGONs are …

invalid

invalid

valid

3.1 – PostgreSQL Installation

- Windows Installer
- PostgreSQL 8.4.4
- Automatically installs itself as a Windows Service

- PgAdmin III
- Psql Interactive SQL Shell

- PostgreSQL 8.4.4

3.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

- Directories created during installation:
- \bin - Executables
- \include - Include files for compilation
- \lib - DLL shared library files
- \share - Extensions

3.1 – PostgreSQL Installation

- Tools included with the install:
- PgAdmin III
- psql Command Line

2.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

3.1 – PostgreSQL Installation

- Under Linux
http://www.postgresql.org/docs/8.4/static/installation.html

- ./configure
- gmake
- su
- gmake install
- adduser postgres
- mkdir /usr/local/pgsql/data
- chown postgres /usr/local/pgsql/data
- su - postgres
- /usr/local/pgsql/bin/initdb -D /usr/local/pgsql/data
- /usr/local/pgsql/bin/postgres -D /usr/local/pgsql/data >logfile 2>&1 &
- /usr/local/pgsql/bin/createdb dbname
- /usr/local/pgsql/bin/psql dbname

3.2 – PostGIS Installation

- As of PostgreSQL 8.3, an Application Stack Builder is used to obtain the latest installers for desired modules.
- Demonstrate Installation …

3.2 – PostGIS Installation

3.2 – PostGIS Installation

3.2 – PostGIS Installation

3.2 – PostGIS Installation

3.2 – PostGIS Installation

3.2 – PostGIS Installation

- The Stack Builder has finished downloading the PostGIS installer.
- Demonstrate PostGIS Installation …

3.2 – PostGIS Installation

3.2 – PostGIS Installation

3.2 – PostGIS Installation

3.2 – PostGIS Installation

3.2 – PostGIS Installation

3.2 – PostGIS Installation

- Under Linux
http://postgis.refractions.net/docs/ch02.html

tar xvfz postgis-1.5.1.tar.gz

cd postgis-1.5.1

./configure

make

make install

createdb yourdatabase

createlang plpgsql yourdatabase

psql -d yourdatabase -f postgis.sql

psql -d yourdatabase -f postgis_comments.sql

psql -d yourdatabase -f spatial_ref_sys.sq

3.3 – Environment Setup

- In order for PostgreSQL to be accessible from the command prompt
- Create a PGHOME environment system variable
- Add %PGHOME%\bin to the system path

3.3 – Environment Setup

3.3 – Environment Setup

3.3 – Environment Setup

3.4 – Create Spatial Database

- From the command line
createdb -T template_postgis dbname

3.4 – Create Spatial Database

- Launch pgAdminIII

3.4 – Create Spatial Database

- Create a new database
- Select “template_postgis” as the template

- Verify the installation
- Check for the spatial system tables
- spatial_ref_sys – stores unique coordinate references
- geometry_columns – stores metadata for a spatial column

- Check for the spatial system tables

3.4 – Create Spatial Database

3.4 – Create Spatial Database

3.5 –Spatially-Enable an Existing DB

- From the command line
psql -f "%PGHOME%\share\contrib\postgis-1.5\postgis.sql" dbname

psql -f "%PGHOME%\share\contrib\postgis-1.5\spatial_ref_sys.sql" dbname

3.5 –Spatially-Enable an Existing DB

- Using a User Interface
- Connect to your existing database
- Open the “Execute Arbitrary SQL Queries” window
- Load/run the PostGIS extension (postgis.sql)
- Load/run the PostGIS spatial reference systems (spatial_ref_sys.sql)

4.1 - Simple Spatial SQL

- “Manually” create geometries
CREATETABLE points

(pt geometry, name varchar);

INSERTINTO points VALUES

('POINT(0 0)', 'Origin');

INSERTINTO points VALUES

('POINT(5 0)', 'X Axis');

INSERTINTO points VALUES

('POINT(0 5)', 'Y Axis');

SELECT name, ST_AsText(pt),

ST_Distance(pt, 'POINT(5 5)')

FROM points;

4.1 - Simple Spatial SQL

name | st_astext | st_distance

--------+------------+------------------

Origin | POINT(0 0) | 7.07106781186548

X Axis | POINT(5 0) | 5

Y Axis | POINT(0 5) | 5

(3 rows)

4.1.1 - PostGIS Reference

- Function References:
http://postgis.refractions.net/docs/ch07.html

“Appendix A: Data Dictionaries”

- (included in ‘PostGIS Workshop.doc’)
Or online PostGIS Documentation

http://postgis.refractions.net/docs

- (included in ‘PostGIS Workshop.doc’)

4.1.2 - Exercises

- Time for some exercises!

4.2 - Loading Shape Files

- Shape File (Misnomer! 3+ Files!)
- .shp = geometry
- .dbf = attributes (string, number, date)
- .shx = utility index

- PostGIS/PostgreSQL Table
- Columns can be geometry
- Columns can be attributes

- One Shape File = One PostGIS Table

4.2 - Loading Shape Files

- shp2pgsql [opts] shapefile tablename
- shp2pgsql –i –s 3005 –D bc_pubs.shp bc_pubs > bc_data.sql

- Read in .shp file
- Write out .sql file
- Load .sql file into PostgreSQL
- using psql
- using PgAdmin

4.2 - Loading Shape Files

- As mentioned in the install process, to make life easier, include %PGHOME%\bin\ in your system path:
- Run cmd.exe
- > cd %WORKSHOP%\data\
- > dir *.shp
- > shp2pgsql -s 3005 bc_pubs.shp public.bc_pubs > bc_data.sql

4.2 - Loading Shape Files

- > notepad bc_data.sql

4.2 - Loading Shape Files

- Load the data via psql
- > psql –d postgis –U postgres –f bc_data.sql

4.2 - Command Line Options

4.2 - Exercises

- Time for some exercises!

bc_pubs100s of points

bc_roads100s of thousands of lines

bc_hospitals10s of points

bc_municipality100s of polygons

bc_habitat_areas100s of polygons

4.3.1 – Indexes and Query Plans

- Databases are fancy engines for speeding up random access to large chunks of data.
- Query plans are the rules used by databases to convert a piece of SQL into a strategy for reading the data.

4.3.1 – Indexes and Query Plans

- In PostgreSQL, these plans can be viewed by prefacing your query with “EXPLAIN”.
- To view the actual run time along side the estimated query plan, preface with “EXPLAIN ANALYZE”.

4.3.1 – Indexes and Query Plans

- In PgAdminIII, view a simple filtered query by pressing the “Explain” button

4.3.1 – Indexes and Query Plans

CREATEINDEX bc_roads_name_idx ON bc_roads (name);

ANALYZE bc_roads;

4.3.2 – When Query Plans Go Bad

- As previously mentioned, a database attempts to minimize disk accesses.
- Indexes can lower disk accesses for queries that only return a few rows, but can actually increase accesses for queries that return a large number of rows.

4.3.2 – When Query Plans Go Bad

- The database builds “plans” based on statistics about data distribution sampled from the tables
- Always tries to be “selective”, to pull the fewest number of records necessary to move on to the next step

- The database builds bad plans when it has bad statistics
- Make sure your database has up-to-date statistics by running ANALYZE

4.3.2 – When Query Plans Go Bad

- It’s important to tune PostgreSQL’s configuration parameters so it can properly chose the best query plan for any given query.
- IE. Is it faster to perform an index scan from a minimally selective query than it is to perform a sequential scan through a table since most of the data is in memory anyway?

Same bounding box, different selectivity!

4.3.3 - Exercises

- Time for some exercises!

4.4.1 – Creating Spatial Indexes

- PostGIS implements R-Tree indexes on top of the GiST indexing system
- Organizes data into nesting rectangles for fast searching

4.4.1 – Creating Spatial Indexes

- An example index creation command is
CREATEINDEX bc_pubs_gidx ON bc_pubs

USINGGIST (the_geom);

4.4.2 - Exercises

- Time for some exercises!

4.4.3 - Using Spatial Indexes

- Spatial index operator is “&&”
- “Bounding Boxes Interact”

A&&B = TRUE

A&&B = FALSE

4.4.3 - Using Spatial Indexes

- Bounding Boxes are not enough!
- Two pass system required
- Use bounding boxes to reduce candidates
- Use real topological test to get final answer

A&&B = TRUE_ST_Intersects(A, B) = FALSE

4.4.3 - Using Spatial Indexes

ST_Intersects(A, B)

A && B AND _ST_Intersects(A, B)

4.4.3 - Using Spatial Indexes

A && B

4.4.3 - Using Spatial Indexes

A && B

4.4.3 - Using Spatial Indexes

_ST_Intersects(A,B)

4.4.3 - Using Spatial Indexes

- Index operations (&&) are built into common functions for automatic use, but you can use them separately if you like.
- ST_Intersects(G1,G2)
- G1 && G2 AND _ST_Intersects(G1,G2)

- ST_Contains(G1,G2)
- ST_Within(G1,G2)
- ST_Touches(G1,G2)
- ST_DWithin(G1,G2,D)
- G1 && ST_Expand(G2,D) AND ST_Distance(G1,G2) > D

- ST_Intersects(G1,G2)

4.4.4 - Exercises

- Time for some exercises!

4.5 - Spatial Analysis in SQL

- A surprising number of traditional GIS analysis questions can be answered using a spatial database and SQL.
- GIS analysis is generally about filtering spatial objects with conditions, and summarizing the results – and that is exactly what databases are very good at.

4.5.1 - Exercises

- Time for some exercises!

4.6 – Data Integrity

- PostGIS permits a user to store invalid or non-simple geometries in a spatial database, unlike Oracle, which require validity upon insertion.

4.6 – Data Integrity

- However, most spatial predicates or functions in PostGIS expect the data to conform to the OGC Simple Features for SQL standard
- Polygon rings don’t cross other rings or self-intersect
- Interior rings on polygons can touch exteriors, but only at a point
- Multi-polygons are non-overlapping

4.6 – Data Integrity

- PostGIS offers ST_IsSimple() and ST_IsValid() functions, useful in determining if geometries conform to the specification.
SELECT count(*)

FROM bc_habitat_areas

WHERENOTST_IsValid(the_geom);

4.6.1 - Exercises

- Time for some exercises!

4.7 – Distance Queries

- It’s easy to write inefficient queries for distances in PostGIS because it’s not always obvious how to use the index in a distance query.

4.7 – Distance Queries

Question: How many wildlife habitats are within 20 kilometers of the municipality of Oliver?

First, find out where Oliver is located.

SELECT ST_AsText(ST_Centroid(the_geom))

FROM bc_municipalities

WHERE name = 'Oliver';

POINT(1470065.29710885 484600.794443135)

4.7 – Distance Queries

Then, use that location to sum all habitats within 20km. The “obvious” way does not use any indexes.

SELECT count(*) AS num_habitats

FROM bc_habitat_areas

WHERE

ST_Distance( the_geom,ST_GeomFromText('POINT(1470065 484600)',3005) ) < 20000;

4.7 – Distance Queries

The “optimized” way is to use ST_DWithin()which silently adds an index into the mix.

SELECT count(*) AS num_habitats

FROM bc_habitat_areas

WHERE

ST_DWithin( the_geom,ST_GeomFromText(

'POINT(1470065 484600)',3005),

20000 );

4.7 – Distance Queries

Here’s the definition of ST_DWithin().

CREATEFUNCTION ST_DWithin(geometry, geometry, float8)RETURNS boolean AS' SELECT $1 && ST_Expand($2,$3) AND $2 && ST_Expand($1,$3) AND ST_Distance($1, $2) < $3 'LANGUAGE'SQL'IMMUTABLE;

4.7.1 – Distance Not Buffer

For complete, here’s how you should never answer a distance query.

SELECT count(*) AS num_habitats

FROM bc_habitat_areas

WHERE

ST_Contains(

ST_Buffer(the_geom, 20000),

ST_GeomFromText('POINT(1470065 484600)',3005)

);

5.1 – OGC Specifications

5.1 – OGC Specifications

- SRID(): Integer - Returns the Spatial Reference System ID for this geometric object. This will normally be a foreign key to an index of reference systems stored in either the same or some other datastore.
- Is3D(): Integer - Returns 1 (TRUE) if this geometric object has z coordinate values.

5.1 – OGC Specifications

- Disjoint(geometry):Integer - Returns 1 (TRUE) if this geometric object is “spatially disjoint” from anotherGeometry.
- Within(geometry):Integer - Returns 1 (TRUE) if this geometric object is “spatially within” from anotherGeometry.
- Contains(geometry):Integer - Returns 1 (TRUE) if this geometric object is “spatially contain” from anotherGeometry.

5.1 – OGC Specifications

- Distance(geometry):Double - Returns the shortest distance between any two Points in the two geometric objects as calculated in the spatial reference system of this geometric object. Because the geometries are closed, it is possible to find a point on each geometric object involved, such that the distance between these 2 points is the returned distance between their geometric objects.
- Intersection(geometry): Geometry – Returns the geometric object that represents the Point set intersection of this geometric object with another geometry.

Approach

make pair-wise tests of the intersections between the Interiors, Boundaries, and Exteriors of two geometries and to represent these relationships in an “intersection” matrix

Possible values:

{T, F, *, 0, 1, 2}

Where:

T == {0,1,2}

F == empty set

* == don’t care

0 == dimensional 0 – point

1 == dimensional 1 – line

2 == dimensional 2 - area

Geometry Topology

Boundary

the set of geometries of the next lower dimension

Point

(dim-0)

Polygon

(dim-2)

Line

(dim-1)

Geometry Topology

Interior

the points that are left when the boundary points are removed

Point

(dim-0)

Polygon

(dim-2)

Line

(dim-1)

Geometry Topology

Exterior

consists of points not in the interior and boundary

Point

(dim-0)

Polygon

(dim-2)

Line

(dim-1)

(b)

(a)

2

1

2

1

0

1

2

1

2

ST_Relate(a, b) = ‘212101212’

5.3 - Exercises

- Time for some exercises!

5.4 – Spatial Predicates

Crosses(geometry):Integer

Crosses relation applies to P/L, P/A, L/L and L/A situations.

Case aP, bL or Case aP, bA or Case aL, bA:

a.Crosses(b)

(I(a) ∩ I(b) ≠∅) (I(a) ∩ E(b) ≠∅)

a.Relate(b, ‘T*T******’)

Case aL, bL:

a.Crosses(b)

dim(I(a)∩I(b)) = 0

a.Relate(b, ‘0********’);

5.4 – Spatial Predicates

Overlaps(geometry):Integer

The Overlaps relation is defined for A/A, L/L and P/P situations.

Case aP, bP or Case aA, bA:

a.Overlaps(b)

(I(a) ∩I(b)≠∅) (I(a) ∩E(b)≠∅) (E(a) ∩I(b)≠∅)

a.Relate(b, ‘T*T***T**’)

Case aL, bL:

a.Overlaps(b)

(dim(I(a) ∩I(b) = 1) (I(a) ∩E(b)≠∅) (E(a) ∩I(b)≠∅)

a.Relate(b, ‘1*T***T**’)

6.1 - Spatial Joins

- Standard joins use a common key
SELECT a.var1, b.var2

FROM a, b

WHERE a.id = b.id

- Spatial joins are based on a spatial relationship
SELECT a.var1, b.var2

FROM a, b

WHEREST_Intersects(a.geom, b.geom)

6.1 - Spatial Joins

- The previous distance query can be expressed using a single spatial join.
SELECT count(*) AS num_habitats

FROM

bc_municipalities a,

bc_habitat_areas b

WHERE

ST_DWithin(

a.the_geom,

b.the_geom,

20000)

AND

a.name ILIKE 'Oliver%';

6.1 - Spatial Joins

Question: Find all pubs located within 250 meters of a hospital

6.1 - Spatial Joins

SELECT bc_hospitals.name, bc_pubs.name

FROM bc_hospitals, bc_pubs

WHEREST_DWithin(

bc_hospitals.the_geom, bc_pubs.the_geom, 250);

name | name

--------------------------------+-----------------------------------

British Columbia Cancer Agency | Holiday Inn Vancouver Centre

The Richmond Hospital | Executive Airport Plaza Hotel ...

St. Paul's Hospital | Holiday Inn Hotel & Suites ...

St. Paul's Hospital | Bosnman's Hotel

St. Paul's Hospital | Sheraton Vancouver Wall Centre ...

St. Paul's Hospital | Burrard Motor Inn

St. Paul's Hospital | Century Plaza Hotel and Spa

St. Paul's Hospital | The Fountainhead Pub

Nanaimo Regional General Hos... | Windward Pub

(9 rows)

6.1 - Spatial Joins

Question: Determine the density of roads per municipality.

6.1 - Spatial Joins

SELECT

m.name,

(sum(ST_Length(ST_Intersection(

m.the_geom, r.the_geom))) /

sum(ST_Area(m.the_geom))) AS road_density

FROM

bc_municipalities m,

bc_roads r

WHERE m.the_geom && r.the_geom

GROUPBY m.name

ORDERBY road_density DESC;

name | road_density

----------------------------+-----------------------

Silverton | 0.000246331953966186

Slocan | 0.000131565932828081

Montrose | 0.000118947787979122

New Denver | 0.000114169586372941

Clinton | 9.92704482763833e-005

. . . . . . .

6.2 - Overlays

- Table-on-table overlays are possible with the ST_Intersection() function
- ST_Intersects(a,b) returns BOOLEAN
- ST_Intersection(a,b) returns GEOMETRY

ST_Intersects() = TRUE

ST_Intersection() =

6.2 - Overlays

Question: Create a new table of habitat areas clipped by the Chilliwack municipal boundary.

6.2 - Overlays

CREATETABLE ch_habitat_areas AS

SELECT

ST_Intersection(h.the_geom, m.the_geom)

AS intersection_geom,

ST_Area(h.the_geom) AS va_area,

h.*,

m.name AS municipality_name

FROM

bc_habitat_areas h,

bc_municipalities m

WHEREST_Intersects(h.the_geom, m.the_geom)

AND m.name = 'Chilliwack';

6.3 – Coordinate Projection

- Every geometry in PostGIS has a “spatial referencing identifier” or “SRID” attached to it.
- The SRID indicates the spatial reference system the geometry coordinates are in.

6.3 – Coordinate Projection

- View the SRID of geometries using the ST_SRID() function
SELECTST_SRID(the_geom)

FROM bc_roads LIMIT 1;

3005

- What’s 3005?
SELECT srtext

FROM spatial_ref_sys

WHERE srid = 3005;

PROJCS["NAD83 / BC Albers",…

- Ah, it’s “BC Albers”

6.3 – Coordinate Projection

- What’s 3005 again?
SELECT proj4text

FROM spatial_ref_sys

WHERE srid = 3005;

"+proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs“

- PROJ4 is coordinate re-projection library used by PostGIS

6.3 – Coordinate Projection

- Coordinate Re-projection is done using the ST_Transform() function
SELECT ST_AsEWKT(the_geom)

FROM bc_roads

LIMIT 1;

"SRID=3005;MULTILINESTRING((

1004687.04355194 594291.053764096,

1004729.74799931 594258.821943696,

1004808.0184134 594223.285878035,

1004864.93630072 594204.422638658,

1004900.50302171 594200.005856311

))"

6.3 – Coordinate Projection

- Coordinate Re-projection is done using the ST_Transform() function
SELECT ST_AsEWKT(ST_Transform(the_geom,4326))

FROM bc_roads

LIMIT 1;

"SRID=4326;MULTILINESTRING((

-125.9341 50.3640700000001,

-125.9335 50.36378,

-125.9324 50.36346,

-125.9316 50.36329,

-125.9311 50.36325

))"

6.4 – Extracting Spatial Data

- One of the utility applications that accompany a typical PostGIS installation is the pgsql2shp command-line utility program that can extract a table, view, or any custom query that contains a geometry column into a shapefile.
USAGE:

pgsql2shp [<options>] <database> [<schema>.]<table>

pgsql2shp [<options>] <database> <query>

6.5 – Clustering

- Clustering is a concept in PostgreSQL where one physically orders the rows in a table based on some index
- CLUSTER [VERBOSE]<tablename> [USING indexname]
- CLUSTER [VERBOSE]

6.5 – Clustering

6.5 – Clustering

7 – PostgreSQL Optimization

- PostgreSQL ships with conservative settings
- Uses very little memory
- Runs on very limited hardware

- Database use cases play an important role in configuring certain parameters
- A read-only webserver with 100s of connections vs. a production system with a few users but high write activity.

7 – PostgreSQL Optimization

- Disk access is slow, so higher performance can be gained by using more memory to cache data!
- shared_buffers
- effective_cache_size

7 – PostgreSQL Optimization

- Sorting is faster in memory
- Increase work_mem

- Disk clean-up is faster with more memory
- Increase maintenance_work_mem

- Allocated per connection
- Also
- Increase checkpoint_segments
- Decrease random_page_cost

8.1 – Basic Exercises

1. What is the perimeter of the municipality of Vancouver?

SELECTST_Perimeter(the_geom)

FROM bc_municipalities

WHERE name = 'Vancouver';

57321.7782018085

8.1 – Basic Exercises

2. How many points were used to digitize the road named “Cartmell Rd”?

SELECTST_NPoints(the_geom)

FROM bc_roads

WHERE name = 'Cartmell Rd';

28

8.1 – Basic Exercises

3. How many holes or interior rings does the habitat area with hab_id = 891 (“Marbled Murrelet”) have?

SELECTST_NumInteriorRings(the_geom)

FROM bc_habitat_areas

WHERE hab_id = 891;

5

8.1 – Basic Exercises

4. What is the length in kilometers of all roads named ‘Main St’?

SELECTsum(ST_Length(the_geom)) / 1000

AS total_length

FROM bc_roads

WHERE name = 'Main St';

44.8489926296202

8.1 – Basic Exercises

5. What is the total area of all municipalities in hectares?

SELECTsum(ST_Area(the_geom)) / 10000

AS total_length

FROM bc_municipalities;

1186482.22780945

8.1 – Basic Exercises

6. What is the average road length (treat every entry in the bc_roads table as a separate road)?

SELECTavg(ST_Length(the_geom))

FROM bc_roads;

434.230644730819

8.2 – Intermediate Exercises

1. How many roads are completely within the municipality of 'Hope‘?

SELECT count(*)

FROM bc_roads r,

bc_municipalities m

WHEREST_Within(r.the_geom, m.the_geom)

AND m.name = 'Hope';

280

8.2 – Intermediate Exercises

2. What is the latitude of the most southerly hospital in the province?

(Hint: consider using a lat/long projection)

SELECT min(ST_Y(

ST_Transform(the_geom, 4269)

)) AS lat

FROM bc_hospitals;

48.4657953714625

8.2 – Intermediate Exercises

3. The residents in the municipality of “Sicamous” should have a watchful eye out for what native wildlife animal?

SELECT b.name

FROM bc_municipalities a,

bc_habitat_areas b

WHERE a.name = 'Sicamous'

ANDST_Intersects(a.the_geom, b.the_geom);

Grizzly Bear

8.2 – Intermediate Exercises

4a. According to the datasets, is there a hospital in “Sicamous”?

SELECT count(*)

FROM bc_municipalities a,

bc_hospitals b

WHERE a.name = 'Sicamous'

ANDST_Contains(a.the_geom, b.the_geom);

0

8.2 – Intermediate Exercises

4b. Is there a pub?

SELECT count(*)

FROM bc_municipalities a,

bc_pubs b

WHERE a.name = 'Sicamous'

ANDST_Contains(a.the_geom, b.the_geom);

1

8.2 – Intermediate Exercises

5. Extract the bc_municipalities polygonal table to a shapefile, except add a new attribute called “perimeter” populated accordingly.

C:\> pgsql2shp -f municipalities.shp -u postgres postgis "SELECT *, ST_Perimeter (the_geom) AS perimeter FROM bc_municipalities"

8.3 – Advanced Exercises

1. What is the name of the closest hospital from Sicamous and how far away is it?

SELECT a.name, b.name,

ST_Distance(a.the_geom, b.the_geom)

FROM bc_municipalities a,

bc_hospitals b

WHERE a.name = 'Sicamous'

ORDER BY ST_Distance ASC

LIMIT 1;

name | name | st_distance

---------+-------------------------------+------------

Sicamous | Shuswap Lake General Hospital | 22943.88318

3.3 – Advanced Exercises

2. What are the names of the pubs located within 100 meters from “Granville St”, located in downtown Vancouver?

SELECTDISTINCT c.name

FROM bc_municipalities a,

bc_roads b,

bc_pubs c

WHERE a.name = 'Vancouver'

AND b.name = 'Granville St'

AND ST_Contains(a.the_geom, b.the_geom)

AND ST_DWithin(b.the_geom, c.the_geom, 100)

ORDER BY c.name;

8.3 – Advanced Exercises

2. What are the names of the pubs located within 100 meters from “Granville St”, located in downtown Vancouver?

name

----------------------

Cecil Hotel

Dufferin Hotel

Fraser Arms Hotel

Howard Johnson Hotel

St. Helen's Hotel

St. Regis Hotel

The Lennox Pub

Yale Hotel

(8 rows)

8.3 – Advanced Exercises

3. What is the total length (in kilometers) of roads named “Douglas St” in Victoria?

(Hint: sum the intersected length)

SELECT sum(ST_Length(

ST_Intersection(r.the_geom, m.the_geom)

)) / 1000 AS sum

FROM bc_roads r,

bc_municipalities m

WHERE r.name = 'Douglas St'

AND m.name = 'Victoria'

ANDST_Intersects(r.the_geom, m.the_geom);

4.47269806771667

8.3 – Advanced Exercises

4. What road would you travel on in BC if you wanted to sight-see BC’s native “Spotted Owl”?

(Hint: try to find the name of the road with the longest intersecting length)

SELECT a.name, sum(ST_Length(

ST_Intersection(a.the_geom, b.the_geom)))

FROM bc_roads a, bc_habitat_areas b

WHERE b.name = 'Spotted Owl'

ANDST_Intersects(a.the_geom, b.the_geom)

GROUPBY a.name

ORDERBY sum DESC;

8.3 – Advanced Exercises

4. What road would you travel on in BC if you wanted to sight-see BC’s native “Spotted Owl”?

(Hint: try to find the name of the road with the longest intersecting length)

name | sum

-----------------------+------------------

Pemberton-Lillooet Rd | 25798.0268482952

Chilliwack Lake Rd | 4937.85768862766

Highway 1 | 4401.56275978918

W Paulsen Rd | 1411.2595495739

E Paulsen Rd | 1226.83009075469

Z005 | 277.194003115705

...

8.3 – Advanced Exercises

5. In the municipality of Kelowna, are the road segments named “Hillaby Ave” disjoint or continuous?

(Hint: consider merging the lines together)

SELECTST_NumGeometries(

ST_LineMerge(ST_Union(the_geom)))

FROM bc_roads

WHERE name = 'Hillaby Ave';

st_numgeometries

------------------

2

Thank You …

- Kevin Neufeldkneufeld@refractions.net
- PostGIShttp://www.postgis.org