Archives de catégorie : Non classé

QGIS : Sélectionner les objets à une certaine distance d’un autre

Un collègue bien embarrassé nous a posé la question suivante : comment sélectionner tous les points d’une couche A qui se trouvent à moins de 4km d’une autoroute.

Le postulat de départ est que nous travaillons sur des fichiers shp (donc on oublie tout de suite la requête sql et l’opérateur st_dwithin() )

Nous avons tous répondu « tu crées un tampon de 4km de rayon autour de l’autoroute puis tu fais une recherche par localisation de tous les objets de la couche de points qui intersectent le tampon créé« .

Il nous a répondu « ok ça je sais le faire mais la personne qui m’a posé cette question m’a tendu un piège pour me dire que c’était fou de devoir passer par un tampon, alors que d’autre outils font ça en un clic…« .
L’objectif sous-jacent et l’effet attendu auprés de l’assistance -non « géomaticienne »- semble être de vouloir déprécier QGIS.

Considérons donc que cette personne n’est pas fan de SQL et qu’elle affectionne le format shp.

Tout d’abord, quel problème y-a-t-il à devoir créer un tampon ?
Aucun… Si les géomaticiens en avaient assez de créer des fichiers à tout va, tous travailleraient en SQL dans une base de données spatiales. Par ailleurs QGIS permet de générer le tampon en mémoire, sans passer par l’écriture d’un fichier sur le disque. Le souci n’est donc pas là.

Peut-être réside-t-il dans le fait de ne pas disposer d’un clique-bouton pour faire le travail et de devoir passer par deux étapes ?
Une des qualités de tout bon géomaticien n’est-elle pas de savoir résoudre un problème qui se pose à lui, avec les outils dont il dispose. Si l’absence de bouton pour faire le travail est un obstacle rédibitoire, changeons tout de suite de métier.

Mais qu’à cela ne tienne. Il me faut un clique-bouton donc je le crée !

Modeleur graphique et boite à outils « traitements »

Le modèle de traitement qui résulte de cette opération est accessible ici en téléchargement.

L’interface produite est la suivante :

Sélection par expression

Une autre possibilité offerte part QGIS réside dans l’utilisation d’expressions pour la sélection.

Je peux donc sélectionner dans une couche de points, tous les objets dont la géometrie est à une certaine distance d’un objet d’une autre couche (ici une ligne).

L’expression utilisée est la suivante :

distance( $geometry , geometry( get_feature( ‘a9′,’nom’,’A9′) ) ) > 4000

Présentation de la solution GeoODK (remplacée depuis par ODK Collect) et de sa mise en oeuvre dans les CEN au « réseau BDD » du CNRS

Présentation réalisée et enregistrée dans le cadre du séminaire « Système d’information embarqué, cahier/carnet de terrain et de laboratoire électronique : quelles interactions avec les bases de données ? »

le Mercredi 05 octobre 2016 à Paris – Jussieu (amphi Charpak).

Cliquez sur l’image ci-dessous pour accéder à la vidéo.

Présnetation d'ODK / Geo ODK

Atelier UICN : les TIC au service de la conservation de la biodiversité

Les démarches présentées sont généralisables et répétables.

La première présentation  concerne des cas d’utilisation d’outils de « reporting » qui croisent connaissance naturaliste et données foncière, pour informer un acquéreur ou un vendeur sur les enjeux connus sur les parcelles concernées ou plus classiquement pour informer un propriétaire de l’intérêt patrimonial de sa propriété :

http://si.cenlr.org/sites/si.cenlr.org/files/UICN/qgis_reporting_carto_chemin.svg

La seconde concerne la possibilité de mettre en place d’outils de saisie de données adaptés au terrain :

  • facilité de création d’un formulaire
  • gain de temps (pas de saisie au bureau)
  • fiabilité (perte d’une session en cas de bug vs. perte d’un carnet de terrain qui couvrent une période beaucoup plus longue)
  • réactivité du système d’information : connaissance immédiate de l’enjeu par les utilisateurs du Si
 
De quoi tester les formulaire et installer la solution : http://si.cenlr.org/geoodk_sicen_mobile
 
Ces deux présentations s’inscrivent dans le cadre plus général du travail des CEN en général que du cas précis du Languedoc-Roussillon.

 

Photographies anciennes de l’IGN – consolidation des PVA dans une table PostGIS

Consolidation par région des geojson générés par cquest dans une table postgis (lambert 93)

-> https://github.com/cquest/photos-aeriennes-ign/tree/master/data/pva
-> http://georezo.net/forum/viewtopic.php?pid=279360#p279360

Les dumps contiennent une table PostGIS ayant la structure suivante :

CREATE TABLE pva_nom_region ( gid integer, idcliche text, mission text, numcli text, idta text, date text, res text, support text, type text, surface double precision, lon double precision, lat double precision, orientation text, url text, largeur integer, hauteur integer, taille integer, geometrie geometry(Polygon,2154) );

Libre à vous de créer la clé primaire et les index :

ALTER TABLE pva_nom_region ADD CONSTRAINT pva_pk PRIMARY KEY(gid); CREATE INDEX pva_geometrie_idx ON pva_nom_region USING gist (geometrie);

La fonction à utiliser dans l’action « ouvrir fichier » pour télécharger ces images depuis QGis (merci David ;-)):

‘http://wxs.ign.fr/clé_d_acces/jp2/DEMAT.PVA/’|| »mission »||’/’|| « url »

Le dump pour PostGIS (France entière toutes années sauf 1941) est ici : pva.sql_.tar.

La couche au format kpkg (France entière toutes années sauf 1941) est ici : pva.gpkg_.tar


Comparer deux versions d’une observation dans la base de données

logo postgresql

 

Cette fonction permet la comparaison deux versions d’une observation.

Elle exploite la table saisie.suivi_saisie_observation avec les fonctions windows : http://docs.postgresql.fr/9.2/functions-window.html

 

La table saisie.suivi_saisie_observation est créée sur le même modèle que la table audit de la documentation plpgsql :

  • http://docs.postgresql.fr/9.1/plpgsql-trigger.html#plpgsql-trigger-audit-example

Chaque version de ligne de la table saisie.saisie_observation y est datée (date_operation). La fenêtre est une partition réalisée selon id_obs et ordonnée par date_operation. Les fonctions window sont expliquées ici dans la doc de postgresql :

Des tests sont à ajouter, notamment pour vérifier que les versions de lignes existent.

Et on peut aller plus loin dans le comportement et afficher les valeur avant/aprés des champs modifiés.

Faire un atlas carto avec QGIS, sans géométrie dans la couche de couverture

Suite aux discussions des journées utilisateurs de QGIS des 10 et 11 décembre. J’ai testé l’utilisation d’une table atributaire comme couche de couverture. Et ça fonctionne !

Je ne sais pas si c’était le cas sur les versions précédentes mais ça fonctionne sur la 2.12. Ma table ne contient que deux colonnes cd_ref et lb_nom.

On utilisera les règles de symbologie pour n’afficher que les données correspondant au cd_ref courant :

 « cd_ref »  = attribute(@atlas_feature ,’cd_ref’)

QGIS rocks !

QGIS comme outil de reporting

Voici la présentation faite à l’occasion du séminaire utilisateur de QGIS des 10 et 11 décembre dernier.

Elle a été réalisée avec Inkscape et SOZI.  Utilisez les boutons de la souris ou les flèches de votre clavier pour avancer/reculer. La molette ou les signes +/- pour zoomer/dézoomer.

http://si.cenlr.org/sites/www.cenlr.org/files/users/webmestre/qgis_et_reporting/qgis_reporting_carto_chemin.svg

Le réseau des CEN mécène du séminaire francophone 2015 des utilisateurs de QGIS

Les Cen sont de gros utilisateurs de logiciels libres, notamment dans leurs SIG avec Qgis et PostgreSQL/PostGIS.

Ces outils nous facilite grandement la gestion et l’exploitation de nos données pour l’exercice de nos missions de conservation des especes naturels.

Cette première initiative du réseau est donc un premier, « petit » mais juste, retour vers la communauté.

SiCen Mobile : utilisation de formulaires ODK pour alimenter notre base de données d’observations

Cette page est le résultat d’un travail mené en commun, en mars 2015 par 4 géomaticiens et informaticiens des CEN Rhône-Alpes (Rémy Clément, Guillaume Costes et Laurent Poulin) et Languedoc-Roussillon (Mathieu Bossaert)

Elle a été actualisée le 17 mai 2018.

OpenDataKit est une suite d’outils libres dédiée à la collecte de données sur terminaux mobiles Androïd.

D’une relative simplicité de mise en oeuvre, la solution permet facilement de décrire et créer des formulaires correspondant à nos besoins. Une fois les données récupérées, il est simple de les intégrer à notre base de données en place.

Nous allons donc passer en revue l’installation des outils de la suite, la définition du formulaire avec XLSForm, et la ventilation des données récoltées dans notre base de données « métier », SiCen.

Présentation générale

ODK est un générateur de formulaires Open Source pour Android. Il permet de collecter des données en mode déconnecté. Les données sont envoyées quand une connexion est disponible, ou par upload de fichiers.

Les formulaires sont créés de manière simple, en utilisant un outil dédié (ODKBuild) ou en les décrivant dans un fichier excel avec le standard XLSForm

Tous les types de données sont disponibles et disposent de « widgets » adaptés : dates, textes, nombres, booléens, geo. Tous les médias que peut créer votre appareil androïd peuvent être attachés à l’observation eux aussi : son, vidéo, photo.

Il est possible d’interroger de longs référentiels (ex. TAXREF), fournis en csv avec le formulaire.

Les outils de la suite

Build

Application en ligne qui permet la création de formulaires grâce à une interface : http://build.opendatakit.org/

XLSForm

« XLSForm est une norme de formulaires créée pour aider à simplifier la création de formulaires dans Excel ». C’est dans un tableur que nous allons décrire de manière simple la strcuture et la logique du formulaire.

Collect

C’est l’application Androïd à proprement parler : https://opendatakit.org/use/collect/

Installée depuis le dépot des applications de google, elle va se connecter au serveur « Aggregate », récupérer et proposer au téléchargement la liste des formulaires disponibles. Puis envoyer à « Aggregate » les données collectées.

Aggregate

C’est le serveur, le chef d’ochestre. Il fournit les formulaires, et consolide et diffuse les données reçues : https://opendatakit.org/use/aggregate/

Installé au sein d’un réseau interne, Aggregate peut-être associé à MySQL ou PostgreSQL. Dans notre cas, il sera asoocié à PostGIS et stockera ses données dans le schéma « odk« de notre base de données « SiCen« .

Mise en œuvre au sein de l’intranet

Installation d’aggregate

Elle peut se faire simplement en utilisant une machine virtuelle, diffusée à prix libre depuis cet été : https://gum.co/odk-aggregate-vm

Dans notre cas, nous allons procéder au déploiement d’une application java sur un serveur Tomcat.

Nous considérerons par la suite que tomcat est installé :  https://opendatakit.org/use/aggregate/tomcat-install/

 

Étapes

On télécharge « l’installeur », qui n’en est pas un : https://opendatakit.org/downloads/

Il s’agit en fait d’un exécutable qui va générer l’archive .war et le script SQL de création de la base de données (base de données, utilisateur odk et schema), conformément aux paramètres renseignés. Les tables sont créées au lancement du war et à l’ajout de nouveaux formulaires

On exécute sur le serveur de base de données les commandes SQL générées. Elles vont permettre de créer les tables et autres objets de base de données nécessaires à Aggregate pour fonctionner.

Puis on déploie l’applicatiopn depuis l’interface de tomcat-manager ou en déplaçant le war dans le dossier webapps de tomcat

Mise en œuvre du formulaire

Nous avions le projet de réaliser une application mobile complète dédiée à l’outil WEB SiCen, permettant de saisir nos observations sur un terminal Android pour les retrouver dans l’interface web de SICEN, sans intervention de l’observateur, autre que d’envoyer les données du formulaire au serveur.

Besoin simple

Il nous faut collecter des données d’observations d’espèces, localisées. Elles sont collectées dans le cadre d’une étude et selon un protocole particuliers

-> une étude, un protocole, des localités sur lesquelles on observe des espèces

Un formulaire GeoODK

→ conçu avec XLSForm

Le fichier excel décrivant le formulaire ainsi que les ressources csv nécessaires à son fonctionnement sont présents dans l’archive ci-jointe : demo_aten.zip
Trucs et astuces

    Pour éviter la demande d’ajout de groupe, lors des « repeat » : Supprimer dans le xml le jr template il est présent au début du xml dans la liste des balises champs du début
    Possibilité de cumuler la fonction quick + search dans la colonne « appearance » du fichier excel des formulaires avec la fonction quick search(…)
    Le widget date avec calendrier n’est pas adapté : mettre dans la colonne « appearance » no-calendar
    Génération des noms de formulaire automatiquement en rajoutant une colonne « instance_name » dans l’onglet « settings » du fichier excel. En utilisant un champ calculé.

Des référentiels en csv générés depuis la base de données

ODK permet désormais d’associer au formulaire de grosses listes de références dans des fichiers csv. Elles seront diffusées à ODKCollect avec le formulaire. Ce dernier les transformera en base de données locale sqlite.

La fonction search() permet de filtrer les entrées proposées, à celle contenant la chaîne de caractèressaisie par l’utilisateur. Nous utiliserons cette possibilité pour la gestion des listes de choix relatives aux :

  • taxons observés
  • observateurs et structures
  • études et protocoles

Voici une démonstration du résultat

 

Création du fond de carte

L’application permet non seulement de créer tout type d’objet mais également d’embarquer sur le terminal des fonds cartographiques en local (pas de réseau nécessaire).

Les fonds doivent être au format mbtiles. Ils peuvent être généré avec TileMill, comme décrit dans la doc de GeoODK : http://geoodk.com/mbtiles_howto.php, ou d’autres outils.

Dans notre cas, nous utiliserons MOBAC parceque nous l’utilisons déjà. Par contre, le mbtile produit par MOBAC nécessitera une petite mabnipulation :
Ouvrir le fichier avec un sqlite manager et exécuter la requête suivante :

CREATE VIEW images as SELECT ROWID, « tile_data », tile_row as « tile_id » FROM « tiles » ORDER BY ROWID

Le fichier MBTiles correspondant aux fonds cartographiques doit ensuite être placé dans un sous dossier du dossier OfflineLayers : OfflineLayersSous_Dossier_contenant_Mbtiles.

Ventilation des données dans la base «métier»

Création d’une vue qui met les données au format propre à notre table de destination (saisie.saisi_obsevration)

Voilà donc nos données envoyées au serveur Aggregate, directement visibles et modifiables dans SiCen !

Conclusion / Bilan

Améliorations récentes

→ utilisation de gros référentiels + widget cartographique

Facilité de mise en œuvre de la solution

→ appli Android + déploiement WAR / Machine virtuelle

Souplesse / facilité de création de formulaires de saisie

→ Par des collègues non géomaticien / montée en compétence rapide

Intégration aisée au SI en place dans la structure

→ En utilisant les outils standards de notre base de données (vues et triggers)

  En lien avec cet article

Validation automatique de donnees

Chaque donnée intégrée à la base de données de l’atlas doit être examinée (validée / invalidée). Afin de faciliter le travail de validation, les fonctions présentées ici, permettent de passer chaque donnée saisie au crible des connaissances actuelles sur l’espèce, issues de la base de données.

Chaque taxon a tout d’abord été « caractérisé » selon les connaissances actuelles mobilisables dans la base de données (données validées).

Pour chaque espèce ont donc été calculées les références suivantes :

  • liste des communes
  • liste des entités et des ensembles paysagers
  • lises des semaines et des décades d’observation
  • listes des observateurs ayant déjà une donnée validée pour ce taxon
  • la répartition altitudinale (plages de 100 m d’altitude, correspondant à la division entière de l’altitude par 100)

Pour chaque nouvelle donnée, ou à chaque modification d’une donnée, les valeurs saisies vont être confrontés à la grille précédemment calculée. C’est un trigger qui déclenche cette confrontation de la donnée saisie aux valeurs de référence pour le taxon

La fonction de comparaison qui est appelée par le trigger est la suivante :

Un « score » est affiché dans le champ « décision validation » afin de permettre aux validateurs de filtrer les données selon ce score.

Dés qu’une donnée est validée, la grille est recalculée pour le taxon concerné. Si une donnée anciennement validée est invalidée, les valeurs de référence pour le taxon sont aussi recalculées.

Au niveau de la base de données, c’est encore un trigger qui déclenche la mise à jour de la ligne de la grille relative au taxon mentionné dans la donnée.

Il appelle cette fonction, qui va supprimer la ligne de référence pour le taxon concerné et la réinsérer en tenant compte des nouvelles données validées :

Intégrer les statuts de protection des espèces animales et végétales

L’INPNpublie à nouveau la liste des espèces faisant l’objet de statuts de protection, ainsi que la liste de ces statuts :
http://inpn.mnhn.fr/telechargement/referentielEspece/reglementation

Après avoir téléchargé les deux fichiers csv, on les intègre au schéma inpn de notre base de données sous les noms inpn.protection_especes et inpn.protection_especes_types

Afin de pouvoir exploiter pleinement cette donnée, nous allons créer une matrice contenant une ligne par taxon, et une colonne par statut. Si l’espèce est concernée par le statut, l’article est mentionné.

Selon la version de postgresql dont on dispose (9.3+ ou inférieure à 9.3), nous allons créer une vue matérialisée ou une table.
La vue matérialisée permettra une mise à jour future de la matrice des protections sans avoir à supprimer, recréer et réindexer cette matrice.

Exporter les données d’observation au Format DEE

Pour le congrès de Dunkerque d’octobre 2015 et le travail d’analyse de la contribution des CEN à la connaissance naturaliste, nous avons entrepris de faire remonter les deonnées du réesau au format « DEE » du SINP.

Ce travail nous a permis de consolider prés de 2000000 de données.

Voici la rquête SQL qui nous a peris cette mise en forme depuis la base sicen.

SELECT row_number() OVER (ORDER BY (lpad(saisie_observation.id_obs::text, 7, '0'::text))) AS gid,
	 CASE saisie_observation.determination::text
		 WHEN 'Vu'::text THEN 'te'::text
		 WHEN 'Entendu'::text THEN 'te'::text
		 WHEN 'Indice de présence'::text THEN 'te'::text
		 WHEN 'Cadavre'::text THEN 'te'::text
		 WHEN 'Capture'::text THEN 'te'::text
		 WHEN 'Collection'::text THEN 'Co'::text
		 WHEN 'Littérature'::text THEN 'Li'::text
		 ELSE 'te'::text
	 END AS statutsource,
	 CASE saisie_observation.determination::text
		 WHEN 'Littérature'::text THEN 'à préciser'::text
		 ELSE NULL::text
	 END AS referencebiblio,
	 lpad(lot_donnee.id_lot::text, 4, '0'::text) AS jddid,
	 'SICEN-LR:00-175'::text AS jddcode,
	 lpad(saisie_observation.id_obs::text, 7, '0'::text) AS identifiantorigine,
	 NULL::text AS identifiantpermanent,
	 CASE
		 WHEN lot_donnee.libelle::text ~~* 'bénévolat%'::text THEN 'Pr'::text
		 WHEN saisie_observation.id_etude = ANY (ARRAY[39,
													   51]) THEN 'NSP'::text
		 ELSE 'Ac'::text
	 END AS dspublique,
	 NULL::text AS codeidcnpdispositif,
	 'CEN L-R'::text AS organismestandard,
	 CASE
		 WHEN lower(saisie_observation.type_effectif) = 'absence'::text THEN 'No'::text
		 ELSE 'Pr'::text
	 END AS statutobservation,
	 btrim(concat(saisie_observation.nom_vern, ' / ', saisie_observation.nom_complet), ' / '::text) AS nomcite,
	 saisie_observation.cd_nom::integer AS cdnom,
	 taxref.cd_ref::integer AS cdref,
	 'non'::text AS sensible,
	 COALESCE(saisie_observation.effectif, saisie_observation.effectif_min, 1::bigint)::integer AS denombrementmin,
	 COALESCE(saisie_observation.effectif, saisie_observation.effectif_max)::integer AS denombrementmax,
	 CASE
		 WHEN saisie_observation.type_effectif ~~* ANY (ARRAY['abondance%'::text,
															  'classe%'::text]) THEN 'Es'::text
		 ELSE 'Co'::text
	 END AS typedenombrement,
	 CASE
		 WHEN saisie_observation.type_effectif IS NOT NULL THEN 'In'::text
		 ELSE NULL::text
	 END AS objetdenombrement,
	 CASE
		 WHEN saisie_observation.observateur = '20'::text THEN 'NSP'::text
		 ELSE md.liste_nom_auteur(saisie_observation.observateur, ', '::text)
	 END AS identiteobservateur,
	 CASE md.liste_nom_structure(saisie_observation.structure, ', '::text)
		 WHEN 'Pas de structure'::text THEN 'indépendant'::text
		 ELSE md.liste_nom_structure(saisie_observation.structure, ', '::text)
	 END AS organismeobservateur,
	 'CEN L-R'::text AS organismegestionnairedonnees,
	 CASE
		 WHEN saisie_observation.observateur <> '20'::text THEN md.liste_nom_auteur(saisie_observation.observateur, ', '::text)
		 ELSE NULL::text
	 END AS determinateur,
	 btrim(concat(validateur.nom, ' ', validateur.prenom)) AS validateur,
	 saisie_observation.remarque_obs AS commentaire,
	 COALESCE(saisie_observation.date_obs, saisie_observation.date_debut_obs) AS datedebut,
	 COALESCE(saisie_observation.date_obs, saisie_observation.date_fin_obs) AS datefin,
	 NULL::TIME WITHOUT TIME ZONE AS heuredebut,
	 NULL::TIME WITHOUT TIME ZONE AS heurefin,
	 NULL::date AS datedeterminationobs,
	 saisie_observation.elevation::numeric AS altitudemin,
	 saisie_observation.elevation::numeric AS altitudemax,
	 NULL::numeric AS profondeurmin,
	 NULL::numeric AS profondeurmax,
	 NULL::text AS codehabitat,
	 NULL::text AS refhabitat,
	 st_asgml(saisie_observation.geometrie) AS geometrie,
	 CASE saisie_observation."precision"
		 WHEN 'GPS'::saisie.enum_precision THEN 10
		 WHEN '0 à 10m'::saisie.enum_precision THEN 10
		 WHEN '10 à 100m'::saisie.enum_precision THEN 100
		 WHEN '100 à 500m'::saisie.enum_precision THEN 500
		 ELSE 1000
	 END AS "precision",
	 CASE
		 WHEN st_geometrytype(saisie_observation.geometrie) ~~* '%POINT%'::text THEN 'St'::text
		 WHEN st_geometrytype(saisie_observation.geometrie) ~~* '%POLYGON%'::text
			  OR st_geometrytype(saisie_observation.geometrie) ~~* '%LINE%'::text THEN 'In'::text
		 ELSE 'NSP'::text
	 END AS natureobjetgeo,
	 COALESCE(saisie_observation.code_insee, commune.code_insee::text) AS codecommune,
	 commune.nom AS nomcommune,
	 sites_cen_inpn_2014.id_mnhn::text AS codeen,
	 NULL::text AS typeen,
	 NULL::text AS codemaille,
	 st_x(st_centroid(saisie_observation.geometrie)) AS st_x,
	 st_y(st_centroid(saisie_observation.geometrie)) AS st_y
FROM saisie.saisie_observation
JOIN inpn.taxref_v8 taxref USING (cd_nom)
JOIN ign_bd_topo.commune ON st_intersects(commune.geometrie, saisie_observation.geometrie)
LEFT JOIN referentiels_divers.sites_cen_inpn_2014 ON st_intersects(sites_cen_inpn_2014.geometrie, saisie_observation.geometrie)
LEFT JOIN md.lot_donnee ON saisie_observation.id_etude = lot_donnee.id_etude
AND saisie_observation.id_protocole = lot_donnee.id_protocole
AND ((CASE
          WHEN st_geometrytype(saisie_observation.geometrie) ~~* '%Point%'::text THEN 'point'::text
          WHEN st_geometrytype(saisie_observation.geometrie) ~~* '%LineString%'::text THEN 'ligne'::text
          WHEN st_geometrytype(saisie_observation.geometrie) ~~* '%Polygon%'::text THEN 'perimetre'::text
          ELSE NULL::text
      END || '_'::text) || CASE saisie_observation.regne
                               WHEN 'Plantae'::text THEN 'espece'::text
                               WHEN 'Animalia'::text THEN 'espece'::text
                               WHEN 'Fungi'::text THEN 'espece'::text
                               WHEN 'Habitat'::text THEN 'habitat'::text
                               ELSE NULL::text
                           END) = lot_donnee.type_donnee::text
LEFT JOIN md.personne validateur ON saisie_observation.validateur = validateur.id_personne
WHERE taxref.cd_nom::text ~ '^[\d+]'::text
  AND (saisie_observation.regne = ANY (ARRAY['Animalia'::text,'Plantae'::text,'Fungi'::text]))
  AND (md.liste_nom_structure(saisie_observation.structure, ', '::text) = ANY (ARRAY['CEN LR'::text,'Pas de structure'::text]))
  AND saisie_observation.cd_nom <> '000000'::text;

TP PostGis Raster

Opérateurs, fonctions et requêtes spatiales

Introduction

Le support des données matricielles (raster) est une « nouveauté » de PostGIS 2.

A chaque pixel de la grille est associée une valeur (altitude, température, densité de population, type d’occupation du sol, perméabilité…). Cette valeur est stockée dans une bande.
PostGIS permet d’associer plusieurs bandes à un raster. On peut par exemple disposer d’un raster de températures mensuelles moyennes qui contient une bande par mois de l’année.

PostGIS stocke chaque raster dans une table de la base de données. La vue raster_columns est l’équivalent pour les raster de la vue geometry_columns pour les données vecteur.
Il offre la possibilité de tuiler le raster, c’est à dire de le fractionner en unités plus petites qui permettront d’accélérer les traitement en utilisant des index.
Chacune de ces tuiles un un carré de x pixels de coté. Chaque tuile (tile) orrespond à une ligne de la table. C’est à dire que la table qui contient le raster compte autant de lignes que le raster compte de tuiles.

De nombreuses opérations sont possibles avec les raster. Référez vous à l’aide de PostGIS1 pour avoir une vision complète des possibilités offertes, ainsi qu’à la présentation initiale du format par ses auteurs :

http://trac.osgeo.org/postgis/wiki/WKTRaster/Documentation01

Objectifs et organisation du TP

L’objectif de ce TP est de découvrir le type de données « raster » supporté par PostGIS depuis la version 2.0. Comme pour les données vectorielles utilisées au cours de deux précédents TP, ce TP  va vous initier aux opérateurs, fonctions et requêtes spécifiques à ce type de données. Chacune des questions posées doit être traitée de la façon suivante :

  • Si cela n’est pas donné, recherche du ou des opérateurs et fonctions permettant de répondre à la question dans le cours ou dans la documentation de PostGIS.
  • Construction de la requête et répondre à la question

Pour une liste des opérateurs et fonctions spatiales ainsi que des exemples d’utilisation, vous pouvez vous référer à cette partie de la documentation PostGIS :

http://postgis.net/docs/reference.html

Une « cheat sheet » relative aux fonctions raster de PostGIS est disponible ici :

http://www.postgis.us/downloads/postgis21_cheatsheet.pdf

Pré-requis

Pour ce TP vous utiliserez la même base de données que pour les TP précédents. Nous utiliserons le client graphique PgAdmin pour l’affichage des tables et des résultats de requêtes.

Références et supports en ligne

Sites traitant de l’utilisation des rasters dans PostGIS

http://postgis.net/docs/manual-2.1/using_raster_dataman.html

Ajout des données Raster à la base de données

Données utilisées

Pour ce TP, nous allons utiliser les données raster du projet WorldClim (http://worldclim.org), disponibles à l’échelle de la planète sur internet :
l’altitude les variables bioclimatiques du projet worldclim :

  • BIO5 = Max Temperature of Warmest Month
  • BIO6 = Min Temperature of Coldest Month

Ces fichiers raster sont lourds (1.4 Go par variable), nous avons donc réalisé une extraction des données du territoire français (métropole et Corse). Vous trouverez ces données sur l’ENT.

En complément, nous allons ajouter à notre base de données la couche vecteur de l’ensemble des communes de France issues de la base de données GEOFLA de l’IGN :

http://professionnels.ign.fr/geofla

Utilisez Quantum Gis (ajout d’une couche raster) pour visualiser l’ensemble de ces couches. Les raster vont apparaître uniformément gris. Pour avoir une idée de la variabilité de la valeur des différents pixels, double cliquez sur chacune des couches. Dans l’onglet Style, choisissez « Bande grise unique » et modifiez la palette de couleur en choisissant «Pseudo-Couleurs »

Ajout de la composante communale de GEOFLA à la table commune

Rendez vous dans le répertoire qui contient votre couche shp de la base de données geofla de l’iGN avec la commande « cd »

Shell> cd /home/…/geofla

Dans ce répertoire utilisez la commande shp2pgsql, vue au cours des séances précédentes pour ajouter ce shp à la base de données. La table de destination s’appellera « communes », la colonne géométrique sera nommée « geometrie », la projection utilisée sera le Lambert 93 (code EPSG 2154) et elle sera indexée.

Attention, l’encodage de cette couche est le LATIN1. La commande shp2pgsql dispose de l’option « W » qui vous permet de spécifier cet encodage.

Shell> shp2pgsql -s 2154 -g geometrie -I -W LATIN1 COMMUNE.SHP communes_france | psql -h gis -U login nom_base_de_données

Cette couche de données est en Lambert 93. Pour faciliter les calculs futurs, qui seront effectués sur des raster en WGS 84, nous allons ajouter une colonne géométrique en WGS 84 (nommée the_geom) à notre table des communes, la remplir et l’indexer. Vous pouvez utiliser PgAdmin3 ou psql pour écrire et exécuter ces requêtes.

psql> ALTER TABLE communes_france  ADD column the_geom geometry(‘MULTIPOLYGON’,4326);

psql> UPDATE communes_france  SET the_geom=st_transform(geometrie, 4326);
Mise à jour de la colonne « the_geom » de la table commune avec la géométrie contenue dan sla colonne « geometrie » reprojetée en WGS 84.

psql> CREATE INDEX communes_geom_gist ON communes_france USING gist (the_geom);
Création de l’index sur cette nouvelle colonne

Ajout des rasters à la base de données

Comme pour les données vecteur, PostGIS propose un chargeur de données en ligne de commande qui s’appelle raster2pgsql.

Consultez la documentation en ligne de PostGIS ou utilisez la commande suivante pour savoir quelles sont les options disponibles :

Shell> raster2pgsql – ?

Le système de coordonnées des grilles que nous allons intégrer à la base de données est le WGS 84 (code EPSG : 4326), nous créerons des tuiles de 5 pixels de coté et nous indexerons la table.

Nous intégrerons le raster bio_5_fr.tif dans une table nommée temperature_max, le raster bio_6_fr.tif dans une table nommée temperature_min et le raster d’élévation nommé alt_fr.tif dans une table nommée altitude.

Pour cela, rendez vous dans le répertoire contenant ces données avec la commande cd :

Shell> cd /chemin/vers/mes/donnees/raster
Shell> /chemin/vers/raster2pgsql -s 4326 -I -t 5×5 bio_5_fr.tif temperature_max | psql -h gis -U login nom_base_de_données
Shell> /chemin/vers/raster2pgsql -s 4326 -I -t 5×5 bio_6_fr.tif temperature_min | psql -h gis -U login nom_base_de_données
Shell> /chemin/vers/raster2pgsql -s 4326 -I -t 5×5 alt_fr.tif altitude | psql -h gis -U login nom_base_de_données

Requêtes géographiques avec PostGIS

Requêtes d’information sur les objets géographiques

Quelle est la taille des pixels des raster « temperature_min » et « temperature_max’ ? En quelle unité est-elle exprimée ?

Utilisez les fonctions ST_PixelHeight() et ST_PixelWidth()

De combien de tuiles sont composés chacun des raster ?

Quelle est la température minimale du centroïde de la commune de Montpellier ?

Vous utiliserez les fonctions ST_intersection(rast,the_geom) et ST_Centroïd(geometry)
Est-ce cohérent ? -> http://worldclim.org/faq

Quelle est la température minimale moyenne de la commune de Montpellier ?

Vous utiliserez la fonction ST_intersection(rast, the_geom) et l’opérateur d’agrégation AVG()

Quelle est la commune la plus chaude du département de l’Aude ?

Quelle est la commune la moins élevée du département des Pyrénées Orientales ?

Création et utilisation d’un raster multibande

Nous allons, pour simplifier les calculs utiliser la possibilité d’ajouter des bandes aux pixels du raster.
Commençons par ajouter au raster altitude une bande contenant la température minimale du mois le plus froid :

UPDATE altitude SET rast=ST_AddBand(altitude.rast, min.rast)
FROM temperature_min min
WHERE ST_UpperLeftX(min.rast) = ST_UpperLeftX(altitude.rast)
AND ST_UpperLeftY(min.rast) = ST_UpperLeftY(altitude.rast);

Ajoutons ensuite une bande contenant la température la plus élevée du mois le plus chaud :

UPDATE altitude SET rast=ST_AddBand(altitude.rast, max.rast)
FROM temperature_max max
WHERE ST_UpperLeftX(max.rast) = ST_UpperLeftX(altitude.rast)
AND ST_UpperLeftY(max.rast) = ST_UpperLeftY(altitude.rast);

Actualisons les informations relatives à nos tables contenues dans la tables raster_columns :

SELECT AddRasterConstraints(‘altitude’::name, ‘rast’::name);
SELECT AddRasterConstraints(‘temperature_min’::name, ‘rast’::name);
SELECT AddRasterConstraints(‘temperature_max’::name, ‘rast’::name);

Vérifions le nombre de bande de notre raster d’altitude :

SELECT num_bands FROM raster_columns WHERE r_table_name = ‘altitude’;

Requêtes complexes

Affichez l’altitude moyenne, la moyenne de la température la plus basse du mois le plus froid, la moyenne de la température la plus élevée du mois le plus chaud de la commune de Sanilhac-Sagries

Ajouter aux communes l’information relative à l’altitude de leur centroïde

  •     ajouter un attribut « altitude » de type numeric à la table « communes_france »
  •     mettre à jour la valeur de cette attribut en interrogeant le raster (utilisez la fonction st_value(rast, geometry), le code de la région est 91…
  •     visualisez la couche dans Quantum Gis en utilisez une symbologie mettant en évidence la variation d’altitude

Bonus pour aller plus loin

Extraire du MNT les zones basses (altitude < 50 cm) des communes littorales de la région de Montpellier (PALAVAS,  LATTES, MAUGUIO, PEROLS, LA GRANDE-MOTTE, VILLENEUVE LES MAGUELONE, VIC LA GARDIOLE, FRONTIGNAN)

On extrait tout d’abord la géométrie et la valeur de chacun des pixels de la table altitude qui sont contenus dans les communes littorales

SELECT communes_france.nom_comm, (st_intersection(altitude.rast, communes.the_geom)).val AS val,
(st_intersection(altitude.rast, communes.the_geom)).geom AS geom
FROM communes_france, altitude
WHERE communes_france.nom_comm::text ILIKE ANY(array[‘PALAVAS%’,’LATTES%’,’MAUGUIO%’,’PEROLS%’,’%GRANDE-MOTTE%’,’%MAGUELONE’,’%GARDIOLE’,’FRONTIGNAN’]) AND st_intersects(altitude.rast, communes_france.the_geom)

De cet ensemble de géométries et de valeurs, nous allons extraire et regrouper les géométrie des pixels de moins de 50 m d’altitude

WITH valeurs_tous_pixels AS (
SELECT communes_france.nom_comm, (st_intersection(altitude.rast, communes_france.the_geom)).val  AS val, (st_intersection(altitude.rast, communes_france.the_geom)).geom AS geom
            FROM communes_france, altitude
           WHERE communes_france.nom_comm::text ILIKE ANY(array[‘PALAVAS%’,’LATTE%’,
‘MAUGUIO%’,’PEROLS%’,’%GRANDE-MOTTE %’,’%MAGUELONE’,’%GARDIOLE’,’FRONTIGNAN’])
AND st_intersects(altitude.rast, communes_france.the_geom)
        )
SELECT valeurs_tous_pixels.nom_comm, avg(valeurs_tous_pixels.val) AS moyenne_temp_min, st_union(valeurs_tous_pixels.geom) AS the_geom
FROM valeurs_tous_pixels
WHERE valeurs_tous_pixels.val < 50::double precision
GROUP BY valeurs_tous_pixels.nom_comm;

On crée la vue.

CREATE OR REPLACE VIEW zone_inondable AS …

Solution proposée

-- 3.1.3

SELECT nom_comm,
--(ST_intersection(rast, st_centroid(the_geom))).val,
ST_value(rast,1,st_centroid(the_geom))
  FROM communes_france, temperature_min
  WHERE nom_comm = 'MONTPELLIER'
  AND st_intersects(rast,st_centroid(the_geom))     -- pour rendre la requête plus rapide
  ;

-- 3.1.4

SELECT nom_comm, avg(val)/10 AS temp_min_moyenn
FROM (
SELECT nom_comm,
(ST_intersection(rast, the_geom)).val
  FROM communes_france, temperature_min
  WHERE nom_comm = 'MONTPELLIER'
  AND st_intersects(rast,the_geom)
  ) AS valeur_des_differents_raster         -- Les sous-requêtes doivent être "nommées"
GROUP BY nom_comm;

/* Alternative en utilisant les CTE (Comman table expression)
   http://www.postgresql.org/docs/9.0/static/queries-with.html */

WITH valeur_des_differents_raster AS (
SELECT nom_comm,
(ST_intersection(rast, the_geom)).val
  FROM communes_france, temperature_min
  WHERE nom_comm = 'MONTPELLIER'
  AND st_intersects(rast,the_geom)
  )
SELECT nom_comm, avg(val)/10 AS moyenne_temp_min
FROM valeur_des_differents_raster
GROUP BY nom_comm;

-- 3.1.5
SELECT nom_comm, avg(val)/10 AS moyenne_temp_max
FROM (
SELECT nom_comm,
(ST_intersection(rast, the_geom)).val
  FROM communes_france, temperature_max
  WHERE code_dept = '11'
  AND st_intersects(rast,the_geom)
  ) AS valeur_des_differents_raster       
GROUP BY nom_comm
ORDER BY 2 DESC LIMIT 1;

-- 3.1.6
SELECT nom_comm, avg(val) AS moyenne_altitude
FROM (
SELECT nom_comm,
(ST_intersection(rast, the_geom)).val
  FROM communes_france, altitude
  WHERE code_dept = '34'
  AND st_intersects(rast,the_geom)
  ) AS valeur_des_differents_raster       
GROUP BY nom_comm
ORDER BY 2 ASC LIMIT 1;

-- 3.3.1
SELECT nom_comm, avg(altitude) AS moyenne_altitude, avg(temp_min)/10 AS moyenne_temp_min, avg(temp_max)/10 AS moyenne_temp_max
FROM (
SELECT nom_comm,
(ST_intersection(the_geom, rast, 1)).val AS altitude,
(ST_intersection(the_geom, rast, 2)).val AS temp_min,
(ST_intersection(the_geom, rast, 3)).val AS temp_max
  FROM communes_france, altitude
  WHERE lower(nom_comm) = lower('Sanilhac-Sagries')
  AND st_intersects(rast,the_geom)
  ) AS valeur_des_differents_raster       
GROUP BY nom_comm
ORDER BY 2 ASC LIMIT 1;

-- 3.3.2
ALTER TABLE communes ADD COLUMN altitude numeric;

UPDATE communes SET altitude = ST_VALUE(rast, 1, ST_CENTROID(the_geom))
FROM altitude WHERE ST_INTERSECTS(rast, ST_CENTROID(the_geom));