Advanced Spatial Analysis with Open Source Tools

For Windows and Mac users

• Install PostgreSQL from the EntrepriseDB website
• When prompted to install stackbuilder, accept to install it and choose to install PostGIS under spatial extensions. When prompted, choose to create a spatial database and name it qcbs_workshop. Answer “Yes” to questions about GDAL_DATA, POSTGIS_GDAL_ENABLED_DRIVERS, POSTGIS_ENABLED_OUTDB_RASTERS.

sudo apt-get install postgresql postgis

To create a spatially enabled database

CREATE DATABASE workshop;
\c workshop;
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;

Click here to download the datasets that will be used for the exercises below. Note that all files are in Latitude/Longitude (WGS84 datum).

• obis_toothed_whales.shp → Occurrences of toothed whales (Odontoceti) downloaded from the OBIS website (http://iobis.org/).
• ETOPO1_Bed_g_geotiff_LR.tif → Bedrock based topography downloaded from here: http://www.ngdc.noaa.gov/mgg/global/global.html
• EEZ_IHO_union_v2.shp → Marine and land zones: the union of world country boundaries and Exclusive Economic Zones Boundaries (EEZ) downloaded from here http://www.marineregions.org/downloads.php
• TM_WORLD_BORDERS-0.3.shp → World political borders downloaded from here http://thematicmapping.org/downloads/world_borders.php
• meow_ecos.shp → Marine Ecoregions of the World (MEOW) downloaded from here http://www.marineregions.org/downloads.php

• Introduction to PostGIS
• PostGIS reference

1. Click on Layer…Add Layer…Add PostGIS layer.
2. Click on New to create a new connection
3. Give the connection a name “qcbs workshop”, specify “localhost” as the host, “workshop” as the database, and put your PostgreSQL Username (usually “postgres”).
4. Make sur the DB Manager plugin is activated in the list of plugins.
5. Go to Database… DB Manager. Click on PostGIS and navigate to the spatially enabled database you created earlier.
6. Click on the Import Layer/File icon and choose file called “muni_s.shp”. Click on “Create spatial index” and leave other options unchecked. Choose “muni_s” as the table name. Remember that PostgreSQL table names should never include uppercase letters, spaces, special characters, or accented characters.
7. Repeat these steps to import the “ProtectedAreas-Quebec.shp”, “Hydrography-Surface.shp” and

sum(expression) aggregate to return a sum for a set of records
count(expression) aggregate to return the size of a set of records
ST_GeometryType(geometry) returns the type of the geometry
ST_SRID(geometry) returns the spatial reference identifier number of the geometry
ST_X(point) returns the X ordinate
ST_Y(point) returns the Y ordinate
ST_Length(linestring) returns the length of the linestring
ST_StartPoint(geometry) returns the first coordinate as a point
ST_EndPoint(geometry) returns the last coordinate as a point
ST_Area(geometry) returns the area of the polygons
ST_Perimeter(geometry) returns the length of all the rings
ST_AsSVG(geometry) returns SVG text

Sum the area of all seas

SELECT iho_sea, sum(ST_Area(eez_iho.geom)) as area FROM eez_iho GROUP BY iho_sea ORDER BY area DESC

Extract latitudes and longitudes of each occurrence for “Delpinus delphis”. Note: we have to use a subquery and the ST_Dump function in this case because the whale occurrence data was imported as a MULTIPOINT data and ST_X/ST_Y expects POINT data.

SELECT ST_X(geom), ST_Y(geom) FROM (SELECT (ST_Dump(geom)).geom  FROM obis_whales WHERE tname='Delphinus delphis') a

ST_Contains(geometry A, geometry B): Returns true if and only if no points of B lie in the exterior of A, and at least one point of the interior of B lies in the interior of A.

ST_Crosses(geometry A, geometry B): Returns TRUE if the supplied geometries have some, but not all, interior points in common.

ST_Disjoint(geometry A , geometry B): Returns TRUE if the Geometries do not “spatially intersect” - if they do not share any space together.

ST_Distance(geometry A, geometry B): Returns the 2-dimensional cartesian minimum distance (based on spatial ref) between two geometries in projected units.

ST_DWithin(geometry A, geometry B, radius): Returns true if the geometries are within the specified distance (radius) of one another.

ST_Equals(geometry A, geometry B): Returns true if the given geometries represent the same geometry. Directionality is ignored.

ST_Intersects(geometry A, geometry B): Returns TRUE if the Geometries/Geography “spatially intersect” - (share any portion of space) and FALSE if they don’t (they are Disjoint).

ST_Overlaps(geometry A, geometry B): Returns TRUE if the Geometries share space, are of the same dimension, but are not completely contained by each other.

ST_Touches(geometry A, geometry B): Returns TRUE if the geometries have at least one point in common, but their interiors do not intersect.

ST_Within(geometry A , geometry B): Returns true if the geometry A is completely inside geometry B

Count the occurrences of each species of whales found in Canadian waters.

SELECT a.tname, count(*) as num FROM obis_whales a, eez_iho b WHERE ST_Within(a.geom,b.geom) AND b.country='Canada' GROUP BY a.tname ORDER BY num DESC

SELECT a.* FROM obis_whales a, tm_world b WHERE ST_DWithin(a.geom,b.geom,1) AND b.name='Canada

Doing it with a distance in meters

ALTER TABLE obis_whales ADD column geog GEOGRAPHY(MULTIPOINT,4326)

UPDATE obis_whales SET geog=geography(geom)

ALTER TABLE tm_world ADD column geog GEOGRAPHY(MULTIPOINT,4326)

UPDATE tm_world SET geog=geography(geom)

SELECT a.* FROM obis_whales a, tm_world b WHERE ST_DWithin(a.geog,b.geog,1) AND b.name='Canada

CREATE TABLE ce03_gbif_occ AS SELECT a.*, count(b.geom) as num_occ FROM cer03 a LEFT JOIN gbif_mammals b ON (ST_Within(a.geom,b.geom)) GROUP BY a.id 

CREATE TABLE routes_test AS SELECT (ST_Dump(ST_Intersection(section_vote_31h5_utm.geom,routes_utm.geom))).geom as geom FROM section_vote_31h5_utm, routes_utm WHERE tri_cep='CHATEAUGUAY'

SELECT sum(ST_Length(geom)) FROM (SELECT (ST_Dump(ST_Intersection(section_vote_31h5_utm.geom,routes_utm.geom))).geom as geom FROM section_vote_31h5_utm, routes_utm WHERE tri_cep='CHATEAUGUAY') a

SELECT avg(elevation) FROM elevation_utm a, section_vote_31h5_utm b WHERE ST_Within(a.geom,b.geom) AND tri_cep='OUTREMONT'

raster. A raster layer.
vector. A vector layer.
table. A table.
number. A numerical value. A default value must be provided. For instance, depth=number 2.4.
string. A text string. As in the case of numerical values, a default value must be added. For instance, name=string Victor.
boolean. A boolean value. Add True or False after it to set the default value. For example, verbose=boolean True.
multiple raster. A set of input raster layers.
multiple vector. A set of input vector layers.
field. A field in the attributes table of a vector layer. The name of the layer has to be added after the field tag. For instance, if you have declared a vector input with mylayer=vector, you could use myfield=field mylayer to add a field from that layer as parameter.
folder. A folder.
file. A filename.

Generate bar graph showing the area of each sea located within a distance of 0.5 degrees from the coast of a chosen country

1 - Open the processing toolbox, click on R Script…Tools…Create New R Script.

2 - Copy and paste the following code into the R Script and save it as Countries_Zones_Plot

##QCBS Workshop=group
##Layer=vector
##showplots
data2=tapply(Layer$new_area,Layer$IHO_Sea,mean)
par(mar=c(10,4,4,2))
barplot(data2, las=2)

3 - Open the Graphical Modeler

4 - Specify “Country Seas” as the Model name and “QCBS Workshop” as the group.

5 - Add a Vector Layer to the Model, name it “Country” and specify Polygon as shape type.

6 - Add a Strong to the Model, name it “Country Name”.

7 - Add “Extract by attributes” to the model (Under QGIS geoalgorithms). Specify Country as the Input Layer, = as the operator, NAME as the Selection Attribute and select Country_Name as the Value.

8 - Add a Fixed Distance Buffer to the model. Specify Output from Extract by Attributes as the Input, and set a Distance of 0.5. Choose to dissolve the output.

9- Add another Vector Layer to the model. Name it Zones. Set Polygon as the shape type.

10 - Add Intersection to the model. Intersect between the output of the Buffer and the Zones layer.

11 - Add Field Calculator to the model. Specify new_area as the Result Field Name and $area as the Formula.

12 - Now Add your R script from Step 2 (Countries_Zones_Plot) to the model. The input is the output from 11 (Field Calculator), Call the output Country_Seas.

Your model should now look like this:

13 - Now your model is ready to be run! First Save the Model and call it Country_Seas. Then close the Graphical Modeler and double click on your Model in the Toolbox under Models…. Specify TM_World_Borders as the Countries and EEZ_IHO_union_v2 as the Zones. Specify the name of a Country (ex. Canada).

14 - When the algorithm finishes running, you should see a graph showing the area of each sea along the coast of the chosen country.