Author Archives: ymr

Le map-matching : Raccrocher des points à un axe avec Postgis/Qgis

Dans le cadre de mon boulot j’ai eu à accrocher des points  de relevés GPS sur une axe d’un référentiel routier,  composé de plusieurs polylignes . Il existe des méthodes avec la base de données Spatialite mais le plus simple reste l’utilisation de Potgis, la cartouche spatiale de Postgres.

Cette technique de raccrochage « intelligent » (au plus proche du moins) s’appelle le map-matching et est couramment utilisé par les grand acteurs de l’information géographique (IGN, Google, TéléAtlas, etc..).

Typiquement les erreurs de mesures des GPS font que parfois, au grès de la réception plus où moins bonne du signal, elles s’éloignent de la cible. Or, les erreurs de relevés n’étant jamais partout pareilles,  les points de sont pas à la même distance de l’endroit réel de la mesure.

Dans notre exemple les mesures du relevé doivent être accrochées à l’axe du référentiel routier pour un traitement linéaire ultérieur.   On utilise donc le map-matching pour projeter la mesure du point sur la ligne du référentiel le plus proche partant du postulat que ce doit être, à quelques grosses erreurs de réception près,  majoritairement le cas.

Principe du map-matching

1- Importation des données sous Postgis

Une fois la base de données installée (voir documentation) Qgis permet d’importer rapidement toute donnée dans Postgis grâce aux outils d’importation du gestionnaire de BD.

Import d’une couche Qgis dans Postgis

 

Attention veillez à bien créer un index spatial et une clef primaire lors de l’importation, surtout si le volume de donnée est important.

Importation dans Postgis avec création d’une clef primaire et d’un index spatial

2- Accrochage/projection des points à la ligne la plus proche avec les fonctions st_interpolate_point, st_line_locate_point et st_distance :

Pour effectuer notre map-matching, nous allons utiliser trois  fonctions Postgis :

  • la  première est st_line_locate_point ( géométrie du point, géométrie de la ligne) : elle permet de renvoyer la fraction de ligne la plus proche d’un point
  • la deuxième fonction est st_interpolate_point (geométrie de la ligne, fraction de ligne ou le point projeté est) : elle permet de renvoyer les coordonnées d’un point de la fraction de ligne projetée.

On cumule donc les deux fonction comme ceci :

st_interpolate(géométrie de la ligne, st_line_locate_point (géométrie de la ligne, géométrie du point))

Principe des fonctions d’interpolation d’un point sur une ligne sous Postgis

 

Problème :  si on l’utilise tel qu’elle elle, la fonction ne marche qu’avec une ligne donnée. Or nous volons accrocher les points à leurs routes les plus proches respectives.

Du coup la requête se complique car on va cumuler cette fonction de projection avec une fonction qui va réduire la projection aux arcs de routes les plus proches via une troisième fonction Postgis

  • la troisième fonction est st_distance(géométrie à, géométrie b) : elle permet de donner la distance entre deux géométries.

On  restreint ainsi  la recherche de la route à une distance à 15 m par exemple autour des points du relevés GPS.

Pour combiner les deux on va utiliser la commande JOIN ON  :

Voilà pour le principe, nous allons profiter aussi pour embarquer les données de nos points GPS mais aussi les croiser avec les nom des routes les plus proches.

Attention !! Qgis demande une clef primaire pour pouvoir afficher des couches spatiales que nous allons créer automatique avec la fonction  row_number() OVER ()::integer

Attention aussi !!  pour plus de clarté, mettre les noms de base  en entête de base et  les nom de tables en entête de champs

ex:  osm est ma base et arc_rd_only est ma table ça donne osm.arc_rd_only pour appeler cette table

la requête donne donc ceci :

 

3- Résultat dans une vue et export vers d’autres formats via Qgis

Nous allons construire une vue qui jointe les deux tables . Une vue permet de fixer notre requête comme une table et permet ensuite  la conversion vers d’autres formats mais aussi l’affichage sous Qgis.

Pourquoi une vue matérialisée ? parce que c’est plus rapide ! En effet Postgis peut stocker les vues comme une table afin d’accélérer le calcul de la requête.

Si on créer une vue classique, à chaque appel de la vue, le calcul de la requête se refait à chaque fois….

 

Voilà donc la création de la vue matérialisée :

Et voilà le résultat. En rouge les points GPS, en vert leur résultante sur l’axe.

Les quelques doublons présents sont le résultat de notre choix d’un tampon de 15m autour des points GPS. A affiner donc, sachant toutefois qu’un chiffre plus bas étant lui aussi source d’erreurs (on peut ne pas « trouver » la route la plus proche).

 

Résultat de l’accrochage (en vert)

 

 

Tagged

[Qgis] un séparateur de milliers avec l’éditeur de fonction de Qgis

Qgis est un logiciel très avancé mais pas exempt de lacunes. Une de celles-ci est l’impossibilité de représenter les données numériques avec un séparateur de milliers.

En attendant que Qgis le permette par défaut , il est possible de pallier à ce manque grâce aux fonctions pythons embarquées (depuis la 2.8).

Dans la carte suivante, qui représente la population en Asie, la donnée numérique est trop longue pour une lisibilité satisfaisante.

On ne sépare pas d’un coup d’œil les pays qui comptent des milliers de ceux qui compte des millions (voir milliards) d’habitants.

impossible de distinguer facilement l’ordre de grandeur d’une population lorsqu’on n’a pas de séparateur de milliers

Il est possible bien sûr en python, comme dans n’importe quel langage,  de séparer automatiquement les milliers.

La fonction ci dessous le permet :

 

Pour faire simple celle fonction sépare les chiffres en groupe tous les 3 caractères et les ré-assemble par groupe de trois avec un espace entre eux.

Notez que celle double boucle sur une ligne est une liste en intention, qui fait la beauté du Python.

Notez en outre que le paramètre (val) est passé en chaîne de caractère (string)  pour que la fonction marche. Elle se base en effet sur un calcul de longueur de chaîne len() pour établir les positions des espaces.

On pourrait utiliser d’autres façons de faire (utiliser le modulo par exemple).

Bref cette fonction peut être intégrée à Qgis moyennant quelques règles

1-  il faut ajouter les bibliothèques de python de qgis => pyqgis :


 

2- il faut ajouter une ligne @qgsfunction() qui déclare le nombre d’arguments que l’on va passer à la fonction (ici un seul) :


 

3-  Attention !! dans pyqgis les valeurs d’un champ d’une couche sont des listes, tableaux ou tuples même s’il s’agit d’une valeur unique :

d’où le val[0] lors de la transformation en chaîne de caractère

 

4- Résultat – La fonction donne alors :


 

Maintenant ou peut l’utiliser comme une fonction native dans l’éditeur de fonction de qgis et ensuite partout où l’on a accès à l’utilisation d’expressions.

l’éditeur de fonction de qgis permet de charger des fonction python dans qgis

 

Dans notre cas il s’agit des étiquettes :

utilisation de la fonction « custom » dans l’éditeur d’étiquettes

 

et voilà le travail !  on peut désormais discrétiser rapidement les état en fonction de leur population.

résultat concluant. notre fonction marche bien

Pas beau le « tuning  » Qgis ? 🙂

Tagged , ,

[Qgis – Fonctions avancées du composeur d’impresssion ] Utiliser l’éditeur de fonction python pour afficher des sommes

Le composeur d’impression de Qgis a longtemps été son point faible, notamment par rapport à ceux des produits d’Esri.

Depuis quelques versions cependant il s’étoffe, notamment sur le volet Atlas (que l’on verra prochainement).

Pour autant certaines fonctions ne sont pas (encore) implémentées, comme par exemple la possibilité d’afficher la somme d’une donnée, d’une couche dans le composeur d’impression.

Heureusement depuis la version 2.8, un éditeur de fonction python permet de charger ses propres fonctions dans l’éditeur d’expression. Grasse à cette possibilité nous allons pouvoir remédier à ce manque.

1- Une somme vous avez dit une somme ?

Mon problème est le suivant : j’ai une couche des communes d’un  département avec la population 2016.

J’ai dupliqué la couche avec un tri pour avoir celles qui dépassent 10 000 hab. en rouge.

Dans le composeur, je voudrait afficher la somme de la population totale .

explication_compos

Dans le composeur d’impression, impossible de sommer une donnée

 

Or si on peut afficher dans le composeur une table d’attribut avec les champs et leurs valeurs, impossible d’en faire la somme.

2- un éditeur de fonction python pour les gouverner tous

Heureusement  il y a dans la calculatrice de champs  un onglet « éditeur de fonctions » qui permet d’insérer du code python.

l’éditeur est assez simple :

  • un onglet à gauche qui les recensent tous
  • un onglet à droite avec un éditeur python intégré

sommechamps1

Nous allons donc charger un fonction appelé sommechamps  qui va effectuer l’opération. Cette fonction est donc en python et utilise la bibliothèque python de qgis : Pyqgis.

notez que l’on passe en paramètre de fonction le nom de la couche (layername) et ne nom du champ (fieldname).

Une fois la fonction chargée on l’appelle simplement dans l’onglet expression avec les paramètres voulu comme un autre fonction qgis native. Ici il s’agit de la couche Toutes les Communes et du champ popul_2016.

sommechamps2

Et voilà la résultat si on l’intègre dans une étiquette du composeur de Qgis.

sommechamps3

On peut utiliser cette fonction pour faire la somme d’une table d’attributs comme celle de nos 5 communes de plus de 10 000 habitants par exemple.

sommechamps4

Bon à savoir si vous voulez renommer, supprimer des fonctions, ou encore les éditer avec un autre éditeur python :  celles-ci sont stockées dans le répertoire .qgis2 (sous linux) ou dans le local data/qgis2 sous ouiouin

Notez la compilation en bytecode python (.pyc) qui s’effectue automatiquement quand on charge le fichier de code python (.py) dans l’éditeur de fonction de Qgis.

sommechamps_url

3- Dans les étiquettes aussi…

Notez enfin que ces fonctions sont utilisable partout où la calculatrice de champs sévit, comme dans les étiquettes par exemple.

Ici on va utiliser la fonction  somme totale des habitants du département – celle du champ population de la commune :

sommechamp-etiq1

et voilu le résultat !

un peu abstrait comme exemple mais il peut servir !! 🙂

sommechamp-etiq2

 

 

 

Changer l’encodage par lots de fichiers .shp de l’ISO-8859-1 vers l’UTF-8

Je sais pas vous mais moi, à chaque nouvelle livraison de la BDTopo, je me demande pourquoi avoir encodé en ISO8859-1 (ou cp1252), bref en encodage windows !

En attendant que l’IGN ai l’illumination et encode ses fichiers shapes en UTF8 comme tout bon chrétien, j’ai fait un script shell qui les converti par lot, amen !

il faut :

  • Installer iconv et ogr/gdal au préalable
  • puis placer le script à la racine de l’arborescence avant de l’exécuter

je laisse les windowsiens traduire dans leur sabir MS-DOS (si c’est possible d’ailleurs….)

Puisque j’y suis pour tous ceux qui ne veullent plus de données attributaires accentuées, j’ai un script python qui fait la même chose : convertir en utf8 puis enlever les accents

Pour info  j’ai utilisé les modules python de gdal/ogr.

Peut être l’idée est adaptable avec les modules python de qgis :  pyqgis pour faire un plugin non ? @ suivre

 

 

 

 

[Fonctions avancées de Qgis – 2ème partie] calculs sur la géométrie et formatage d’étiquettes

Nous avions vu il y a quelques temps comment utiliser les fonctions avancées de Qgis afin de croiser des couches entres elles dans l’éditeur de style de Qgis.

Aujourd’hui nous allons voir l’édition avancée des étiquettes.

Qgis permet depuis longtemps, grâce notamment à l’éditeur d’expression, un réglage très fin des étiquettes et les mises-à-jour successives améliorent encore ces réglages.

 

Afficher des résultats de calculs géométriques dans les étiquettes

Nous avons par exemple la couche commune des P.O. de la BDTopo sur laquelle nous voudrions faire apparaître la distance depuis la préfecture, à savoir Perpignan, dans l’étiquette de chaque commune.

 

couche_commune

Description de la table Commune avec sa donnée attributaire

 

Nous allons utiliser l’expression distance de l’éditeur d’expression de qgis qui permet de sortir la distance (en mètre si système de projection de type Lambert ou degrés si wgs84).

Cette fonction compare la géométrie d’une entité avec un autre de manière très simple :

Ici nous allons prendre la géométrie de l’entité courante,  c’est à dire de chaque commune, grace à l’expression $geometry,  que nous allons comparer à celle de Perpignan.

Comme vu dans la première partie,  il est possible de sortir la géométrie d’une autre couche  grâce aux expressionx geometry(entité) et getFeature(couche,champs,valeur du champ).

Ici notre couche se nomme « commune », notre champ de tri sera la code INSEE de la commune « code_insee » et ‘66136’ pour la valeur du code INSEE de la commune de Perpignan (voir illustration plus haut).

Nous allons sortir la valeur arrondie (pas besoin de 10 décimales)  avec l’expression round()  de la distance trouvée en Km.

On ajoute un peu de texte autour !

distance_à-perpignan

Dans l’éditeur d’expression on peut fusionner (concaténer) plusieurs fonctions ou textes avec la commande || note : on peut aussi utiliser le retour à la ligne avec la commande ‘\n’

et voilà !

distance_sans_centroide

Résultat de la fonction de calcul de distance dans les étiquettes avec l’enveloppe de la commune

 

 

Eh !!! mais pourquoi les communes qui touchent Perpignan sont marquées  à 0 Km de distance ? (ex : Pollestres, Canohès, Toulouges, etc.)

En fait l’expression distance() fonctionne sur l’enveloppe de la commune.

Utilisons donc le centroïde pour plus de précision avec l’expression centroid(géométrie)

 

distance&centroide

Fonction de calcul de distance à partir des centroïdes

Le résultat est mieux non ?

 

Résultat de la fonction de calcul de distance dans les étiquettes avec l’enveloppe de la commune

 

Formater une étiquette selon une condition

Ok tout cela est sympa mais il y a  un truc qui me gène dans le résultat précédent.

En effet pour Perpignan,  l’étiquette affiche bêtement « Commune de Perpignan / 0 Km de Perpignan’

Ça fait tâche vous ne croyez pas ?

Réglons le problème avec la mise en place de conditions.

Comme vu aussi dans la première partie, on peut utiliser l’expression CASE qui fonctionne de la sorte :

Pour la condition nous allons utiliser l’expression attribute(entité, champ de l’entité)  qui renvoie l’attribut d’un champ donné mais aussi $currentfeature qui renvoie l’entité en cours.

La fonction suivante compare donc l’entité de la commune de l’étiquette avec celle, fixe, de Perpignan :

en gros attribute(get_Feature(‘communes’,’code_insee’, 66136 ),’nom’) renvoie le nom de la commune de l’entité dont le code INSSE=66136, à savoir => ‘Perpignan’ et le compare au résultat du nom de la commune de l’entité courante ($currentfeature).

Si la condition est remplie, on affiche juste le nom de la commune de Perpignan, sinon on exécute la fonction de distance.

enlever_perpignan

Affichage de la distance si la commune est différente de la commune de Perpignan

 

Resultats_enlever perpignan

La commune de Perpignan n’affiche plus la distance à 0 Km

 

Formatage conditionnel dans Qgis 2.12

La nouvelle version, Qgis 2.12, nom de code  « Lyon », permet en outre d’appliquer des règles directement sur les étiquettes de manière native (sans passer comme vu plus haut par la fonction CASE).

etiquette_condition_qgis_lyon

Le formatage conditionnel dans Qgis 2.12

 

Ici notre filtre est le même que vu précédemment mais sous la forme :

  • soit vrai =>  attribute ($currentfeature,’nom’) = attribute(get_Feature(‘communes’,’code_insee’, 66136 ),’nom’)=1
  • soit faux => attribute ($currentfeature,’nom’) = attribute(get_Feature(‘communes’,’code_insee’, 66136 ),’nom’)=0

On peut ainsi complexifier largement nos opérations sur les étiquettes.

resulat_etiquettes_condition_qgislyon

Le formatage conditionnel des étiquettes permet de multiples combinaisons

 

Le résultat est la même que dans le point précédent.

 

Filtrage d’étiquettes sur une expression régulière

Ok comme le montre l’image ci dessous il nous reste un problème à régler : lorsque la commune commence par un voyelle, il est malhabile d’utiliser le « de » devant.

Par exemple la commune d’Opoul, s’affiche ‘commune de Opoul’, beurk !!

problemes_voyelles

Problème devant les communes qui commencent par une voyelle

 

Pour remédier à cela et afficher un d’ devant les communes qui commencent par une voyelle,  on va utiliser  la fonction regexp_match(chaîne de caractère, expression régulière).

On va tout d’abord sortir la première lettre du nom de commune avec l’expression substr(chaîne de caractère, position de début, nbr de caractères) et supprimer la casse en la mettant d’office en minuscule avec l’expression lower().

Ensuite on va utiliser regexp_match pour vérifier si cette première lettre est bien une voyelle.

Les expressions régulières permettent de trouver un motif particulier (ou patern)  dans une chaîne de caractères.

Je ne m’attarde pas dessus car cela mérite un cours à part entière, mais sachez que dans Qgis les expression régulière sont celles du langage de programmation python.

allez faire un tour ici ou

 

Bref voici la fonction  :

filtre_voyelles

Si la première lettre du nom de commune est une voyelle alors on affiche Commune d’ , sinon Commune de

 

Bon mais le problème est maintenant que je m’aperçoit que certaines communes sont à particules et commencent par « Le » ou « La » comme par exemple Le Barcarès.

dans ces cas, ma fonction échoue, zut !!

problème_nom_composé

Échec de la fonction dans la cas de communes à particules

On va complexifier un peu si vous le voulez bien !

On va ainsi d’abord tester la présence de particules dans un bloc CASE, puis, si ça correspond,  exécuter un autre contrôle avec la focntion if(condition, opération si oui, opération si non)

 

ouch !! ça pique un peu quand même !!

filtre_nom_composé_voyelles

Double condition CASE puis IF dans le ELSE  !!

 

mais le résultat est là !!

 

résultat

 

Pas beau qgis ?

 

 

 

 

 

Créer de beaux fonds cartographiques avec Qgis & des données libres

Mes professeurs de géographie m’ont toujours dit que la lisibilité d’une carte tenait en grande partie au choix du fond de plan.

Les SIG ne dérogent pas à cette règle.

Or combien de fonds hideux ou trop complexes j’ai vu tourner derrières des SIG ?

Mais rien n’est perdu car avec  des outils et des données libres ont peu déjà faire beaucoup de choses !

 

1- Donner du relief à sa cartographie

Un joli fond en relief et d’un coup les informations que l’on va coller dessus vont ressortir.

Pour cela on va utiliser les MNT (Modèle Numérique de Terrain). Comme vous le savez, le MNT est un raster dont la valeur pixel = donnée d’altitude.

On va commencer par les télécharger dans des sites qui les mettent à disposition. Il en existe de toutes tailles et résolutions…

La NASA et le ministère japonais de l’économie mettent à disposition des MNT sur toute la planète via son programme ASTER (Avanced Spaceborn Thermal Emission and Reflexion radiometer)

La NASA et le ministère de l’économie Japonnais mettent à disposition un MNT sur toute la surface du globe

La NASA et le ministère de l’économie Japonnais mettent à disposition un MNT sur toute la surface du globe

 

 

Une fois le MNT téléchargé (le plus souvent il s’agit du format .asc pour ASCII grid), chargez-le dans Qgis.

 

MNT chargé dans Qgis

MNT chargé dans Qgis

 

Comme dis plus haut, le MNT est un raster dont la bande pixel contient l’élévation. Comme toute couche SIG , on peut donc créer une thématique en fonction de l’altitude.

Il en existe de toutes sorte dans les rampes de couleur de Qgis. On choisit le nombre de classe et on choisit ses valeurs de bornes

Mise en place de la thématique en fonction des niveaux d'altitudes

Mise en place de la thématique en fonction des niveaux d’altitudes

Résultat de la thématique

Résultat de la thématique

 

Ça a déjà plus de gueule non ?

Bon maintenant on va construire le relief en donnant un ombrage en fonction de la pente.

De la pente ? kékidi le monsieur ?

Et  oui si on a la valeur  d’altitude, on peut en déduire la pente par analyse de voisinage. Heureusement il existe des outils qui font ça automatiquement !

Personnellement j’utilise GRASS et son module r.shaded.relief mais gdal fait ça aussi.

On peut y choisir l’azimut du soleil, l’angle qui forme par rapport à l’horizon, etc….

Dans tous les cas allez  chercher ces outils d’ombrage dans le boite à outils de qgis (Traitement => boite à outils) .

 

Module r.shaded.relief de GRASS

Module r.shaded.relief de GRASS

 

Le résultat est bien sympa. On peut jouer avec la transparence de la couche d’ombrage pour rendre un effet de relief à notre thématique.

 

Résultat de l'ombrage mis en transparence sur le MNT

Résultat de l’ombrage mis en transparence sur le MNT

 

2- Un  peu de style !! Dessiner le trait de côte et l’océan

On peut aussi désormais marquer le trait de côte et voire même dessiner l’océan.

Pour cela on va faire ressortir les valeur de pixel du MNT qui sont égale à  0, c’est à dire le niveau de la mer.

On va utiliser la calculatrice raster, qui est un bel outil qui permet de faire des opérations sur les rasters et leur(s) bande(s).

 

la calculatrice Raster permet de sortir uniquement les valeurs d'altitude = 0

la calculatrice Raster permet de sortir uniquement les valeurs d’altitude = 0

Le nouveau raster créer fait ressortir en négatif les côtes. Il ne nous rester plus qu’à extraire le contour de même valeur pour avoir le trait de côte.

GDAL et GRASS font ça mais je vous conseille d’utiliser d’utiliser l’outil SAGA qui donne le résultat le plus propre.

 

Résultat du calcul de raster pour les valeurs null. Le trait de côte ressort en négatif.

Résultat du calcul de raster pour les valeurs null. Le trait de côte ressort en négatif.

 

Module SAGA qui permet d'extraire la courbe de niveau zéro

Module SAGA qui permet d’extraire la courbe de niveau zéro

 

 

Et voilà notre trait de côte. On va pouvoir lui imprimer un style désormais.

trait_de _côte

 

 

Nous allons créer un polygone avec l’option lignes vers polygones. Un fois ce polygone créer, nous allons appliquer un dégradé suivant la forme puis un polygone inversé. Ainsi ces deux effets combinés vont donner un style de dégradé vers le large qui colle bien.

 

polygones inversés

Mis en place du style polygone inversé sur un style dégradé suivant la forme

 

style_ocean

Dégradé suivant la forme avec polygones inversés

 

fond_de_carte2

le résultat

 

Et voilà, en quelques minutes nous avons pu créer un fond de carte.

Personnellement je ne met que le relief en blanc (avec transparence > 65%)  et non la rampe de couleur car sinon la carte est trop chargé comme le montre l’image ci-après

 

fond_de_carte

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Suivre l’impact des feux de forêts par imagerie satellite avec le plugin Qgis SCP (Semi-Automatic Classification Plugin) : Exemple de l’incendie de fin mars 2015 en Ukraine à proximité de la centrale nucléaire de Tchernobyl

La télédétection est la technique qui permet d’acquérir de la donnée sur le territoire sans y mettre un pied. Pratique, même si le géographe vous dira qu’il faut ne jamais faire abstraction du terrain !

Bref quand on parle de télédétection, on pense souvent aux images satellites qui permettent de cartographier de zones reculées et/ou en voie de développement, là ou justement la cartographie est lacunaire. Les satellites, tournant à grande vitesse autour du globe, permettent en outre de couvrir de grandes distances et de revenir régulièrement aux même endroit. Ainsi les images satellites permettent aussi un couverture périodique du territoire. La mission Landsat de la NASA par exemple permet de couvrir toute zone sur terre tous les 16 jours.

Ces images servent donc à l’observation directe (encore que les résolutions spatiales ne sont pas souvent fameuses) mais elles enregistrent surtout les fréquences de lumières au delà du spectre visible par l’œil humain.

La méthode est complexe aussi nous n’aborderons que le principe général qui repose sur le réflectance . En effet les capteurs des satellites divisent les longueurs d’ondes renvoyées par la matière présente sur terre en bandes. La possibilité d’avoir la bande infrarouge et proche infrarouge par exemple permet de discrétiser la végétation.

C’est cette méthode de classification qui permet de réaliser des cartes de couverture des sols, comme par exemple le programme européen Corine Land Cover.

 


Nous allons utiliser les deux intérêts de la télédétection par images satellites, la fréquence de prise de vues et le large spectre des capteurs, pour suivre, presque en direct (nous sommes mi mai) et mesurer l’impact de l’incendie géant qui s’est produit en Ukraine le 28 mars dernier.

Le but de ce tutoriel est de mesurer la surface de l’incendie dans les heures et les jours qui suivent sont départ.

Pour ce faire, nous allons utiliser :

– une source de donnée ouverte, à savoir les images satellites Landsat 8,

– le plugin SCP (Semi-Automatic Classification Plugin) pour Qgis.


1- La donnée : images Landsat8

Le satellite Landsat 8 a été lancé en février 2013. Il couvre la Terre tous les 16 jours avec des images de 185 km x 185 km,en 16 bits, comptant 11 bandes spectrales : 9 dans le visible (8 multispectrales de résolution 30 m ; 1 panchromatique à 15m) et 2 thermiques (60 m).

 

Landsat 8
Operational
Land Imager
(OLI)
and
Thermal
Infrared
Sensor
(TIRS)Lancement de la missionle 11 février 2013
Bandes Longueur d’onde
(micromètres)
Résolution
(mètres)
Band 1 – utltra bleue (pour étude côtiere) 0.43 – 0.45 30
Band 2 – Bleue 0.45 – 0.51 30
Band 3 – Vert 0.53 – 0.59 30
Band 4 – Rouge 0.64 – 0.67 30
Band 5 – Near Infrared (NIR) – proche Infrarouge 0.85 – 0.88 30
Band 6 – SWIR 1 – Infrarouge 1.57 – 1.65 30
Band 7 – SWIR 2 – Infrarouge 2.11 – 2.29 30
Band 8 – Panchromatique 0.50 – 0.68 15
Band 9 – Cirrus (détection de nuages cirrus – application météo) 1.36 – 1.38 30
Band 10 – Thermal Infrared (TIRS) 1 – Thermique 10.60 – 11.19 100 * (30)
Band 11 – Thermal Infrared (TIRS) 2- Thermique 11.50 – 12.51 100 * (30)

 

Pour télécharger ces images, on peut aller sur le site de l’IGN américain, l’USGS (US Géological Survey) : http://earthexplorer.usgs.gov/

mais on verra que l’on peut aussi télécharger directement depuis le plugin scp


 

2- Le plugins SCP pour Qgis

Le plug-in est développé par un chercher italien, Luca Congedo. Il est installable simplement depuis qgis. Il a juste besoin des bibliothèques scientifiques de python (scipy, numpy et matplotlib)  pour fonctionner pleinement. Il existe un riche manuel pour ce plugin mais malheureusement pour les anglophiles, il n’existe pas en français.

Le principe de la méthode de classification est simple est repose sur le principe de l’identification d’une petite zone, la ROI pour Région Of Interest, et de son interpolation à tous les pixels similaires du raster.

Un fois installé , le plugin ajoute 3 barres (docks) à qgis :

1- la scp toolbar qui permet de convertir (corrections atmosphériques et de réflectance) , voire directement télécharger, les images Landsat et les rendre utilisable pour la classification. C’est d’ici aussi que se configure l’application

2- la barre ROI  qui permet de sélectionner des zones et leur affecter une longueur d’onde.

3- la barre classification qui permet d’affecter un algorithme de classification définie avec la région d’intérêt et d’étendre la discrétisation à toute la carte.

On va donc utiliser ses outils dans ce même ordre chronologique

 

2.1 Télécharger et convertir les images Landsat8 avec la scp toolbar :

On peut télécharger les images depuis la site de l’USGS ou directement depuis le plugin après avoir sélection les dates et la zone voulue.

L’outil permet donc de télécharger les images Landsat (de la version 4 à 8)  pour une zone donnée (Area coordinate) que l’on définie sur la carte de Qgis,pour une plage de dates données (ici du 29/04 au 16/05 de 2015).

Notez sur l’image ci-après que l’on peut aussi demander une image avec peu de nuages (max cloud cover) et que l’on peut directement les convertir et les charger dans Qgis

télécharger_landsat8

l’outil de téléchargement des images Landsat

 

Le résultat du téléchargement est une suite de rasters, une par bande,  que l’on va ensuite convertir pour les rendre utilisable à la classification. Notez que l’on peut appliquer des méthodes de correction atmosphériques et thermiques assez simplement.

De plus l’outil va créer automatiquement un raster virtuel (.vrt) qui va nous permettre d’éviter de charger toutes les bandes dans Qgis. c’est pas top ? 🙂

conversion_images_landsat-scp

Outil de conversion des images Landsat

 

2.2 Définir la zone de d’intérêt (ROI)

Tout d’abord nous allons devoir localiser l’incendie en lui même.

Facile ! le panache de fumée permet de le localiser de très haut.

localisation_feux

localisation de l’incendie par son panache de fumées

Néanmoins quand il s’agit de repérer les surfaces brûlées,  c’est une autre paire de manche…

Heureusement, nous l’avons vu plus haut, Landsat8 dispose de bandes Infrarouge. Nous allons donc dans les propriétés de qgis affecter aux bandes RGB :

le Proche Infrarouge pour le Rouge, le Rouge pour le vert et le Vert pour le bleu. C’est le mode dit « RGB 543 »

RGB543

Vue du raster en mode RGB 453

 

et voilà le résultat ! on voit maintenant. apparaître en noir les zones brûlées  Elle est pas belle la vie ?

visual_2_ir

Ensuite, après avoir créer un .shp pour l’occasion (nommé roi.shp) afin de garder en mémoire les zone de classification, nous allons prendre un échantillon de la zone brûlée que nous allons définir comme telle.

roi

Étapes de la création de la zone de classification

 

Le plugin va maintenant associer les pixels de la zone ROI définie pour une longueur d’onde donnée. C’est par cette méthode que nous allons pouvoir extrapoler ensuite le résultat à toute la zone.

longueur_d_onde

signature spectrale de la zone définie

 

2.3 – La classification

En fonction de cette signature nous allons pouvoir classifier. Pour cela il faut cliquer sur « add to signature » dans la fenêtre roi.

La classe (ici « terre brulee ») passe ainsi dans la fenêtre de classification

 

add_to_signature

Ajout de la zone à l’onglet de classification

 

Dans le cas d’une classification complète de type LandCover, il faut répéter l’opération plusieurs fois.

Dans notre cas, nous allons procéder à la classification, après avoir prévisualisé nos choix d’algorithme de classification et de réglage de seuil.

Nous avons choisi l’algorithme  « Spectral Angle Mapping » avec un seuil de 10.

 

On obtient un raster avec comme bande la valeur de classification. Un fois transformée en vecteur, nous allons pouvoir trouver l’étendue de l’incendie.

Pour cela nous allons rajouter au vecteur un champ aire de valeur => $area

 

calcul_aire

calculatrice de champ de qgis pour le calcul de l’aire de chaque entité

Voila, on peut sommer toutes les aires des entités ainsi trouvées. On trouve près de 10 000 ha !!

 

3 – Résultats :

Et voilà , garce aux données ouvertes Landsat et au plugin scp de qgis on peut faire du monitoring sur des catastrophes d’ampleur en quasi temps réel.

En l’occurrence c’est une catastrophe qui n’a pas fait beaucoup de bruit médiatique  mais pourtant on voit pourtant bien que l’incendie est très proche de la centrale endommagée !

incendie_ukraine_29_mars_2015

incendie_ukraine_29_mars_2015_infrarouge

 

Tagged , , ,

Analyse de réseau avec GRASS/QGIS [2eme partie] : Analyse de la pertinence de la sectorisation des collèges publics des P.O

Nous avons vu dans la 1ere partie comment utiliser certains des modules d’analyse réseau de GRASS dans Qgis, soit le module de création de réseau v.net et celui de séparation du réseau selon des coûts de distance v.net.iso.

Nous allons maintenant vérifier la pertinence de la sectorisation des collèges (la fameuse carte scolaire) à l’aune des coûts de déplacements, c’est à dire voir dans quelle mesure cette sectorisation respecte le trajet le plus court vers le collège du secteur ou s’il existe un collège plus proche.

 

Note sur les données utilisées pour l’étude :

– les données du réseau routier proviennent de la BDtopo ® de l’IGN (licence gratuite pour la recherche)

– les données de la population scolarisée par communes pour la tranche d’age de 11-15 ans (collèges) proviennent de la population légale de l’INSEE (recensement 2011/2012)

– les données de la sectorisation et des effectif des collèges proviennent du site de l’Académie de Montpellier.

 


Introduction

Pour la suite de notre étude, il faut prendre en compte deux facteurs importants :

1- La sectorisation se fait sur la base communale, comme le montre la carte de la sectorisation ci-dessous, ce qui en soit est une source d’approximation. En effet les tailles de communes n’étant pas égales entre elles, cet échelon du territoire n’est sans doute pas la zone géographique la plus appropriée pour dessiner des aires de proximité autour des Collèges.

2- Ensuite, et c’est un point important, la proximité géographique n’est pas le seul critère, il y a celui du nombre d’enfants scolarisés qui entre en compte dans cette sectorisation comme le montre la carte de la réparation communale du nombre de collégiens

 

sectorisation_colleges

carte_repartition_communale_collegiens

 

1-Les aires d’influence des collèges sur le réseau routier avec l’utilisation du module GRASS : v.net.alloc

Pour trouver ces « aires d’influence » on peut utiliser le module GRASS qui alloue un sous-réseau en fonction du centre le plus proche : v.net.alloc

Il s’agit, une fois le réseau crée avec le module v.net (voir 1ere partie), d’allouer un sous réseau catégorisé par un ou plusieurs centres (donnée ponctuelle) que représente les établissements.

Cette catégorie peut être filtré (grâce à l’identifiant de l’établissement) ou être de la forme 1-402 dans le cadre de la prise de tous les établissements scolaire ,ou encore 1-30 pour le cas des collèges publics uniquement.

v.net.alloc

Création d’un sous-réseau de même distance par rapport aux centres. Remarquez qu’ici les 402 établissements scolaires des P.O sont pris en compte

 

Le résultat permet de mettre en évidence les réseau qui sont les plus proche des collèges, et ce , pour chaque établissement du département

12-v.net.alloc_resulats

Exemple du résultat de v.net.alloc sur les établissements scolaire du département

 

Nous avons obtenus un réseau linéaire d’influence des collèges (voir carte ci-dessous) mais pour plus de lisibilité nous allons le transformer en zone.

 

resultats

Résultats de l’utilisation du module v.net.alloc sur les collèges publics du département uniquement

 

 2- Affiner & présenter les résultats avec Qgis

Nous allons transformer le réseau en raster afin de pouvoir fondre les zones de même distance dans des polygones car c’est toujours plus parlant, notamment à petite échelle.

Pour cela nous commençons par « rastériser » notre réseau en utilisant comme champ des futurs pixels l’identifiant du collège.

 

rasteriser1

La résolution de chaque pixel est fixée à 500 x500 de manière à fusionner un peu plus notre réseau. Et voilà !

vecteur_vers_raster

Ce raster peut maintenant  être interpolé avec d’autres vecteurs mais nous allons juste le retransformer en polygones, toujours avec Qgis.

Pour cela on utilise la fonction de vectorisation raster -> vecteurs :

 

polygoniser

 

Et voilà nous avons nos zones « d’influence » des collèges.

 

polygone_colleges


 

 3- Résultats

Nous allons croiser les données de la sectorisation des collèges avec celle de nos chemins le plus court pour se rendre au collège le plus proche et nous obtenons cette carte qui synthétise la carte scolaire et en modère sa pertinence.

Nous avons mis de côté Perpignan pour plus de lisibilité et car cela mériterait une carte à part.

 

synthese_pertinence_carte_scolaire

 

 

 

 

Tagged , ,

Analyse de réseau avec GRASS/QGIS [1ere partie] : calcul de l’éloignent aux collèges publics des P.O

Dans ce tutoriel nous allons voir comment utiliser les fonctions d’analyse de réseau de GRASS avec Qgis comme interface graphique.

Le sujet d’étude est la proximité aux services publics d’éducation avec notamment la question du coût de déplacement. Cette question est capitale quand on sait que l’heure du réveil peut s’avérer être un facteur prépondérant sur la réussite scolaire.

Introduction

GRASS GIS est LE logiciel SIG opensource historique. Il comporte des modules d’analyse spatiale très avancés, mais le soucis est qu’il est austère.

Son interface graphique et ses fonctions d’édition cartographiques sont en outre mois souple que Qgis, encore que cela est moins vrai avec la nouvelle version stable : la 7.0.

Aussi il existe un plug-in GRASS pour Qgis qui permet de bénéficier du meilleur des deux mondes.

C’est cette fonctionnalité que nous allons utiliser pour notre étude.

Je ne reviens pas sur l’installation et le paramétrage de GRASS pour Qgis qui est bien détaillé dans la documentation de qgis mais il faut retenir deux facteurs importants qui peuvent dérouter ceux qui utilisent GRASS pour la première fois :

1- GRASS possède sa propre base de données, le  dataset, qui est à la croisée des chemins entre le fichier et la base de donnée. Il faut donc importer toute donnée (raster, vecteur ou autre) dans GRASS avant de pouvoir travailler avec.

2- GRASS travaille sur le principe de région, c’est à dire une emprise géographique, qu’il convient de définir avant de commencer.

 


1- Préparation de la donnée

Nous cherchons à analyser les temps de transport autour des établissements scolaires des P.O, en particuliers les collèges publics.

Pour cela nous allons utiliser les données ouvertes du ministère et télécharger la géolocalisation des établissements scolaire sur le site d’opendata de la mission etalab :  data.gouv.fr

Il s’agit d’un fichier délimité (.csv)  qui contient les adresses et le coordonnées des établissements mais aussi des informations supplémentaires (nom, niveau, nature, etc.)

le pdf ci-dessous décrit cette donnée

description_donnee_etablissements_scolaires

 

Nous allons donc commencer par importer cette donnée csv dans qgis grâce à l’outil d’import de csv :

1-import_csv

Import du csv dans qgis. l’import est entièrement paramétrable (encodage, choix du délimiteur, choix des champs x,y, etc.)

 

le résultat donne une carte où les établissements sont représentés par un ponctuel :

2-import_france

Résultat de l’import du .csv

Seulement nous ne voulons travailler que sur le département des P.O, nous allons donc trier cette donnée  grâce à la commande filtre de Qgis.

Nous avons remarqué que la numérotation du ministère comporte le code département au début.

Nous allons donc utiliser l’expression ‘like ‘066%’ qui permet de trier tous les numéro commençant par 066 :

3-tri_sur_po

tri sur la donnée

 

 


2- Import des données dans GRASS

Nous allons désormais importer les points de localisation des établissements scolaires ainsi que le réseau routier de la BDtopo® de l’IGN dans GRASS grâce au module d’import vectorielle v.in.ogr

 

4-import_dans_grass

import de données vectorielles dans GRASS

Le résultat est visible dans l’outil GRASS de qgis qui recense en outre les différents modules :

5-couches_importees

les couches importées et leurs statistiques basiques dans l’outil GRASS pour Qgis

 

 

Un fois importées dans GRASS les couches sont  utilisable avec toutes les fonctions de qgis (thématique, composeur, etc.)

C’est là que réside l’atout majeur de la fusion de ces deux outils.

6-resultat_import

Les couches GRASS dans Qgis (on voit en rouge la région de travail de GRASS)

 


3- Les modules d’analyse réseau de GRASS (v.net.*)

Tous les modules d’analyse réseau de GRASS commencent par v.net (v pour vecteur, net pour réseau). Mais dans tous les cas il faut commencer par affecter un réseau liant les points et le réseau  avec la commande v.net

7.1grass_commande_analyse_reseau

Les modules d’analyse réseau de GRASS

 

3.1 – Création du réseau – module v.net

 Le réseau que l’on va créer va permettre de mettre en liaison les ponctuels représentant les établissements scolaires avec le réseau routier de la BDtopo®

Cette opération est nécessaire avant toute analyse réseau dans GRASS.

Pour cela, on utilise le module v.net . La couche de vecteur principale est le réseau routier, la couche de ponctuels est celle des établissements scolaires.

7-commande_v.net_creation_reseau

création d’un réseau avec GRASS

Nous allons choisir l’option « connecter des points non connectés par insertion de nouvelles lignes » avec un seuil (Treshold) de 500 m

Il s’agit de relier les établissements scolaires au réseau routier, y compris si le point représentant l’établissement est à 500 m du réseau routier.

creation_reseau

le module v.net permet de connecter les points au réseau

 

3.2 – Analyse du réseau selon le coût de déplacement autour des collèges publics du département : module v.net.iso

Nous allons tout d’abord trier les établissements dans qgis pour ne retenir qu les collèges public puis réimporter le tout dans GRASS (et créer un nouveau réseau avec le module v.net)

13-tri_colleges

Tri des collèges publics dans Qgis

 

Une fois le réseau créé, nous allons utiliser le module v.net.iso pour sortir les coûts de déplacements (ici en distance) autour des collèges publics du département.

Mais nous pouvions faire de même avec le temps de déplacement si nous avions la vitesse dans la donnée attributaire du réseau routier.

isolignes

le module v.net.iso de GRASS

 

Ici les paramètre à passer au réseau sont

la catégorie => ici il s’agit de filtrer sur un des 30 collèges ou les prendre tous en compte (1-30)

les coûts (cost)  => ici on choisit les distances que l’on veut étudier (moins de 1km à plus de 30 km dans notre exemple)

4- Résultat 

Et voilà la carte permet de mettre à jour les coûts de déplacements et ainsi la proximité au service public d’éducation.

Mais nous verrons dans d’autres articles d’autres fonctions  toutes aussi puissantes.

@ suivre…

synthese_colleges_cout_distance

 

 

 

Tagged , , ,

[Fonctions avancées de Qgis – 1ère partie] Croiser des couches par comparaison spatiale

Dans tout SIG qui se respecte, et Qgis ne déroge pas à la règle bien sûr,  on peut croiser les couches entre elles.

Par exemple on peut trier la couche Bâti Industriel de la BD TOPO® de l’IGN avec une commune donnée.

Seulement le résultat de ce croisement va produire une troisième couche.

Or il est possible avec les fonctions de Qgis de réaliser ce tri directement depuis la couche voulue.

En l’occurrence si l’on veut afficher uniquement le bâti industriel sur la commune de Perpignan comme illustré ci-après :

Description de la couche Commune de la BD TOPO® de l’IGN

Pour cela on va d’abord récupérer la géométrie de la couche Commune et qui a un CODE_INSEE=66136 avec la commande suivante :

Ensuite on compare cette géométrie avec la géométrie de la couche courante du Bâti industriel => $geometry

Ici on va utiliser la commande within (géométrie a, géométrie b) qui prend les entitées de la géométrie a contenue dans la géométrie b

Si elle le résultat de la commande est égale  à 1 c’est qu’elle est à l’intérieur de la commune de Perpignan :

Nous allons donc mettre cette formule dans la gestion des Styles d’affichage de Qgis (‘Onglet Style’ -> ‘Ensemble de règle’) sur la couche à trier, à savoir la couche bâti industriel.

On obtient uniquement le bâti industriel qui est à l’intérieur de la commune de Perpignan !  La classe non ?

Dans la règle de filtrage de Qgis, le résultat de la fonction avancée permet de faire un tri sur une comparaison spatiale

Attention l’utilisation de cette fonction demande du temps de calcul à chaque déplacements/zoom !!

 


 

Mais on peut aller plus loin et afficher en vert le bâti industriel qui n’est pas sur Perpignan, et en rouge celui qui est sur Perpignan.

il suffit de rajouter une règle où la commande précédente est égale à 0

et voilà la travail !!

 

 


Allons encore plus loin ! oui oui  😉

Je veux afficher en vert l’étiquette du type de bâti hors Perpignan et en rouge celle dans Perpignan comme c’est déjà le cas du polygone du bâti.

Dans l’option couleur de l’étiquette (voir ci-dessous) j’utilise une fonction bien sympa  : la fonction color_rgb

Avec une condition (fonction CASE) on met en rouge (utilisation de l’encodage RGB classique) -> [ color_rgb(255,12,0) ]  lorsque on est sur Perpignan :

 

Et pour la gloire je vais ajouter sur les étiquettes des bâtiments industriels de Perpignan la phrase « sur la commune de Perpignan »

En effet une commande Qgis permet de récupérer le champ d’une couche => c’est la commande attribute (objet, nom du champ à afficher)

Dans mon cas je vais récupérer le champ du nom de la commune (‘NOM’) :

dans l’affichage de la valeur de l’étiquette je vais concaténer le nom de la commune avec la nature du bâtiment :

 

 

Le résultat  :

 

Vous pouvez vous amuser avec ces fonctions qui sont décrite ici

 

Bonne expérimentation !

Tagged , , ,