Atelier du CSBQ sur l'utilisation des logiciels libres pour l'analyse spatiale avancée.

23 février 2017, Université du Québec à Rimouski.

• Lien vers la présentation Prezi
• Présentation Prezi simultanée

Pour utilisateurs Windows et Mac

• Installez PostgreSQL depuis le site web EntrepriseDB
• Quand on vous demande si vous voulez utiliser stackbuilder, acceptez et choisissez d'installer PostGIS sous “spatial extensions”. Quand on vous le demande, choisissez de créer une base de données spatiale et nommez-là spatialdb. Répondez oui aux questions à propos de GDAL_DATA, POSTGIS_GDAL_ENABLED_DRIVERS, POSTGIS_ENABLED_OUTDB_RASTERS.

Pour utilisateurs Linux (Ubuntu/Debian)

sudo apt-get install postgresql postgis

Then type

sudo gedit /etc/postgresql/9.5/main/postgresql.conf

enlevez le # devant #listen_addresses = 'localhost'

redémarrer postgresql

sudo service postgresql restart

Pour créer une base de données spatiale

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

Cliquez ici pour télécharger les jeux de données qui seront utilisés pour l'atelier. Notez que tous les fichiers sont en Latitude/Longitude (datum WGS84).

• obis_toothed_whales.shp → Occurrences des odontocètes (toothed whales) téléchargées depuis le site OBIS(http://iobis.org/).
• world_bedrock_topo.tif → Topographie de la surface terrestre téléchargée ici: http://www.ngdc.noaa.gov/mgg/global/global.html
• EEZ_IHO_union_v2.shp → L'union des zones marines et des “Exclusive Economic Zones Boundaries” (EEZ) téléchargé ici http://www.marineregions.org/downloads.php
• TM_WORLD_BORDERS-0.3.shp → Frontière politiques mondiales, téléchargées ici: http://thematicmapping.org/downloads/world_borders.php
• meow_ecos.shp → Marine Ecoregions of the World (MEOW), téléchargé ici [http://www.marineregions.org/downloads.php]]

• Introduction to PostGIS
• PostGIS reference

1. Cliquez sur Layer…Add Layer…Add PostGIS layer.
2. Cliquez sur New pour créer une nouvelle connexion
3. Donnez un nom à la connexion “atelier_csbq”, spécifiez “localhost” comme hôte, “spatialdb” comme base de données spatiale (ou un autre nom que vous lui avez donné) et mettre votre nom d'usager PostgreSQL (habituellement “postgres”).
4. Assurez-vous que l'extension DB Manager est activée dans la liste des extensions (plugins).
5. Allez à Database… DB Manager. Cliquez sur PostGIS et naviguez vers la base de données spatiales (spatialdb) que vous avez créé plus tôt.
6. Cliquez sur l'icône Import Layer/File et choisissez le fichier “obis_toothed_whales.shp”. Choisissez “Create spatial index” et “Create single-part geometries instead of multi-part” et laissez les autres options inchangées. Nommez la nouvelle table “obis_whales” en minuscules. Souvenez-vous que les noms de tables PostgreSQL ne devraient jamais inclure de lettres majuscules, d'espaces, de caractères spéciaux ou d'accents.
7. Répétez ces étapes pour importer EEZ_IHO_union_V2.shp (nommez là eez_iho), meow_ecos.shp (nommez là meow) et TM_WORLD_BORDERS.shp (nommez là tm_world).

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

Créer une nouvelle couche, seulement avec les bélugas

CREATE TABLE belugas AS SELECT * FROM obis_whales WHERE tname='Delphinapterus leucas'

Superficie de chaque pays, classée en ordre décroissant et calculée en kilomètres carrés.

SELECT name, ST_Area(geography(geom))/1000000 as areakm2 FROM tm_world ORDER BY area DESC

Superficie de chaque mer, classée en ordre décroissant.

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

Extraire les latitudes et longitudes de chaque occurrence de “Delpinus delphis”. Note: nous devons utiliser une sous-requête et la fonction ST_Dump dans ce cas car les occurrences d'odontocètes ont été importées en format MULTIPOINT et les fonctions ST_X/ST_Y exigent des POINT.

SELECT ST_X(geom), ST_Y(geom) FROM obis_whales WHERE tname='Delphinus delphis';

Distance au plus proche voisin pour chaque localisation de béluga

SELECT DISTINCT ON(a.id) a.id, b.id as nn, ST_Distance(a.geom,b.geom) as dist FROM obis_whales a, obis_whales b WHERE a.id <> b.id AND a.tname='Delphinapterus leucas' AND b.tname='Delphinapterus leucas' ORDER BY a.id,dist

EXERCICE 1 Trouvez le pays qui a la plus gros ratio entre la racine carrée de son aire et son périmètre. Quel pays a le plus petit ratio?

Réponse

SELECT name, sqrt(ST_Area(geom))/ST_Perimeter(geom) as ratio FROM tm_world ORDER BY ratio DESC

EXERCICE 2 Extraire les coordonnées X et Y dans le SCR WGS 84/World Mercator du seul individu de “Electra” dans le jeu de données.

Réponse

SELECT ST_X(ST_Transform(geompt,3395)) as X, ST_Y(ST_Transform(geompt,3395)) as Y FROM obis_whales WHERE tname='Electra'

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

Extraire les noms des odontocètes et le nom de la mer dans lesquels ils se trouvent en eaux canadiennes.

SELECT a.tname, b.iho_sea FROM obis_whales a, eez_iho b WHERE ST_Within(a.geom,b.geom) AND b.country='Canada'

Compte les occurrences de chaque espèce d'odontocètes dans les eaux canadiennes.

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

Sélectionne toutes les occurrences d'odontocètes à moins de 1 degré des limites terrestres canadiennes.

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

Pour utiliser une distance en mètres

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

Créer un nouveau polygone avec une zone tampon de 100 km autour de l'Australie.

CREATE TABLE aus_buffer AS
SELECT ST_Buffer(geography(geom),100000)::geometry(Multipolygon,4326) AS geom
FROM tm_world
WHERE name='Australia'

EXERCICE 3 - Calculer le nombre d'occurrence d'odontocète par mois dans les eaux britanniques (table eez_iho). Utiliser la fonction EXTRACT de PostgreSQL pour extraire le mois de la colonne datecollec. Pour convertir une colonne de date textuelle en date, utiliser ::date (datecollec::date).

Réponse

SELECT EXTRACT(MONTH FROM datecollec::date) as month, count(*) as cnt FROM obis_whales a, eez_iho b WHERE ST_Within(a.geom,b.geom) AND b.country='United Kingdom' GROUP BY MONTH ORDER BY month

EXERCICE 4 Dans les pays qui font moins de 20,000 km² (2e+11 m²), lequel a le plus grand nombre d'espèces d'odontocètes à moins de 1 degré de distance de ses frontières?

Réponse

SELECT b.name, count(DISTINCT a.tname) as num_sp FROM obis_whales a, tm_world b WHERE ST_Area(geography(b.geom))< 2e+11 AND ST_DWithin(a.geom,b.geom,1) GROUP BY b.name ORDER BY num_sp DESC

Défi Créez une nouvelle table avec toutes les occurrences d'odontocètes situés dans la région latidudinale polaire(Lat_Zone) au nord de l'équateur. Visualisez ces observations avec la projection North_Pole_Azimuthal_Equidistant dans QGIS.

Réponse

CREATE TABLE whales_polar as SELECT a.* FROM obis_whales a, meow_ecos b WHERE ST_Within(a.geom,b.geom) AND lat_zone='Polar' AND ST_Y(a.geompt)>0

Téléchargez le fichier R contenant le code pour cette partie de l'atelier.

• Package RQGIS
• Package spgrass6
• Package RPostgreSQL
• Package rgdal

• QGIS Training manual. Guide du module Traitements.
• QGIS User manual - Processing framework.

Commandes pour spécifier les fichiers entrants et sortants et les options des scripts Processing

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.

Créez un script Processing/R créant 100 points distribués aléatoirement sur la terre dans R et qui sont importés comme couche dans QGIS.

Réponse

##Atelier CSBQ=group
##points=Output vector
xy=cbind((runif(100)-0.5)*360,(runif(100)-0.5)*180)
data=data.frame(id=1:100)
points <- SpatialPointsDataFrame(coords=xy,data=data)
proj4string(points) <- CRS("+init=epsg:4326")

Créez un diagramme en barre montrant l'aire de chaque mer située à moins de 0.5 degrés de la côte d'un pays choisi.

1 - Ouvrez le “processing toolbox”, cliquez sur R Script…Tools…Create New R Script.

2 - Copiez et collez le code suivant dans le script R et sauvez-le sous Countries_Zones_Plot

##Atelier CSBQ=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 - Ouvrez le “Graphical Modeler”

4 - Specifiez “Mers pays” comme le nom du modèle et “Atelier CSBQ” comme le groupe.

5 - Ajoutez une couche vectorielle au modèle, nommez-là “Pays” et spécifiez polygone comme “shape type”.

6 - Ajoutez une “String” au modèle et nommez-là, “Nom du pays”.

7 - Ajoutez “Extract by attributes” au modèle (sous QGIS geoalgorithms). Spécifiez “Pays” comme “Input Layer”, = comme opérateur, NAME comme “Selection Attribute” et “Nom du pays” comme “Value”.

8 - Ajoutez un “Fixed Distance Buffer” au modèle. Spécifiez “Output from Extract by Attributes” comme “Input”, et mettre la distance à 0.5 (c'est en degrés puisque la carte est en degrés). Choisissez de “dissolve the output”.

9 - Ajoutez un autre vecteur au modèle. Nommez-le “Zones”. Spécifiez “Polygon” comme “shape type”.

10 - Ajoutez “Intersection” au modèle. Spécifiez une intersection entre “output of the Buffer” et la couche “Zones”.

11 - Ajoutez Field Calculator au modèle. Spécifiez new_area comme “Result Field Name” et $area comme Formule. Ceci calculera l'aire des polygones de chaque section de mer dans la zone tampons autour du pays en degrés.

12 - Maintenant, ajoutez le script R créé plus haut à l'étape 2 au modèle. “Input” est la sortie de l'étape 11 (Field Calculator). Appelez le fichier sortant “Pays_Mers”

Votre modèle devrait ressembler à ceci:

13 - Vous êtes maintenant prêt à rouler le modèle! Sauvez d'abord le modèle et appelez le Pays_mers. Ensuite, fermez le “Graphical Modeler” et double-cliques sur votre modèle dans le Toolbox sous Models. Spécifiez TM_World_Borders comme Pays and EEZ_IHO_union_v2 Pour les Zones. Spécifiez le nom d'un pays de votre choix (ex. Canada).

14 - Quand l’algorithme cesse de rouler, vous devriez voir un graphique montrant l'aire de chanque mer le long de la côte du pays choisi.

Créer un modèle graphique dans Processing qui permet d'extraire l'élévation et la pente pour chaque point d'occurrence de deux espèces d'odontocètes choisies et de créer des histogrammes montrant la distribution des pentes et des élévations des deux espèces dans R. Utilisez la fonction r.slope de GRASS et la fonction “Add grid values to points” de SAGA (dans Shapes…points) pour arriver à vos fins.

Aide

##Guillaume analysis=group
##Layer1=vector
##Layer2=vector
##showplots
par(mfrow=c(2,2))
hist(Layer1$slopetif,main=paste('Pente pour',Layer1$tname[1]))
hist(Layer1$worldbedroc,main=paste('Elevation pour',Layer1$tname[1]))
hist(Layer2$slopetif,main=paste('Pente pour',Layer2$tname[1]))
hist(Layer2$worldbedroc, main=paste('Elevation pour',Layer2$tname[1]))