====== Atelier du CSBQ sur l'utilisation des logiciels libres pour l'analyse spatiale avancée. ====== 23 février 2017, Université du Québec à Rimouski. * [[https://prezi.com/gtw8-gwn8gqt/atelier-du-csbq-sur-lutilisation-des-logiciels-libres-pour-lanalyse-spatiale-avancee/|Lien vers la présentation Prezi]] * [[http://prezi.com/gtw8-gwn8gqt/present/?auth_key=txah9hq&follow=GuillaumeLarocque&kw=present-gtw8-gwn8gqt&rc=ref-19631691|Présentation Prezi simultanée]] ===== Installation des logiciels ===== ==== Installation de PostGIS ==== Pour utilisateurs Windows et Mac * Installez PostgreSQL depuis le site web [[http://www.enterprisedb.com/products-services-training/pgdownload|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; ==== Données pour les exercices ==== {{:atelier_specialspatial_fichiers.zip|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]] ===== Utiliser PostGIS avec QGIS ===== * [[http://workshops.boundlessgeo.com/postgis-intro/|Introduction to PostGIS]] * [[http://postgis.net/docs/manual-2.1/reference.html|PostGIS reference]] ==== Importer des données de QGIS dans PostGIS ==== - Cliquez sur Layer...Add Layer...Add PostGIS layer. - Cliquez sur New pour créer une nouvelle connexion - 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"). - Assurez-vous que l'extension DB Manager est activée dans la liste des extensions (plugins). - Allez à Database... DB Manager. Cliquez sur PostGIS et naviguez vers la base de données spatiales (spatialdb) que vous avez créé plus tôt. - 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. - 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). ==== Fonctions de bases pour des requêtes sur des couches uniques ==== **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 | Plus gros: Nauru Plus petit: Canada 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' ++++ ==== Relations spatiales entre couches. ==== **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 [[http://www.postgresql.org/docs/9.3/static/functions-datetime.html|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 km2 (2e+11 m2), 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 ++++ \\ \\ ===== Utilisation de PostGIS, QGIS et GRASS dans R ===== {{::special_spatial_code_r.r|Téléchargez le fichier R contenant le code pour cette partie de l'atelier}}. * Package [[http://cran.r-project.org/web/packages/RQGIS/index.html|RQGIS]] * Package [[http://cran.r-project.org/web/packages/spgrass6/index.html|spgrass6]] * Package [[http://cran.r-project.org/web/packages/RPostgreSQL/index.html|RPostgreSQL]] * Package [[http://cran.r-project.org/web/packages/rgdal/index.html|rgdal]] ===== Utilisation de Processing dans QGIS ===== * [[http://docs.qgis.org/2.14/fr/docs/training_manual/processing/index.html|QGIS Training manual. Guide du module Traitements. ]] * [[http://docs.qgis.org/2.14/en/docs/user_manual/processing/index.html|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.\\ ===== Exercice 5 ===== 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 | Votre fichier de script R devrait ressembler à ceci: ##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") ++++ ===== Exercice 6 ===== 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 "Extracted (attribute)" comme "Input", et mettre la distance à 0.5 (c'est en degrés puisque la carte est en degrés). À "dissolve the output", choisissez 'No'. 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 "output of the Buffer" comme Input et "Zones" comme couche d'intersection. 11 - Ajoutez Field Calculator au modèle. Spécifiez Intersection comme Input. 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: {{::graphical_modeler3_processing.png?direct&200|}} 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-cliquez 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. France). 14 - Quand l’algorithme cesse de rouler, vous devriez voir un graphique montrant l'aire de chaque mer le long de la côte du pays choisi. ===== DÉFI ===== 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 | Votre Script R devrait ressembler à ceci ##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])) Votre modèle graphique devrait ressembler à ceci: {{::temp.png?direct&100|}} ++++