Tag Archives: qgis

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 , ,

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 , , ,

Réaliser un carte d’images aériennes comparées avec QGIS, GDAL et Leaflet

Depuis quelques temps, l’Institut national de l’information géographique et forestière met à disposition gratuitement sur le géoportail de l’IGN des images aériennes de différentes époques.

Le but de ce tutoriel est de créer une cartographie comparée entre les images de la campagne de photos aériennes de 1924 et celles d’aujourd’hui, le tout uniquement à l’aide d’outils libres

 


1- Télécharger les images 

La première chose à faire est de télécharger les images qui nous intéresse sur le Géoportail (voir les étapes 1, 2, et 3 sur l’image ci-après).

Note : Plusieurs années sont disponibles. Il est aussi possible de télécharger de vielles cartes (Cassini et cartes d’État major)

ortho1924_screen1

Choix des séries à télécharger sur le Géoportail

 

ortho1924_screen2

Télécharger les images sélectionnées

 


2- Orienter les images

  1. Premier hic, les images sont au format compressé  Jpeg2000. Ce format est lisible avec gdal via la bibliothèque OpenJpeg mais pas par beaucoup d’autres de logiciels. Sous windows, l’IGN met l’outil IGNMap à disposition pour les lire
  2. Deuxième hic, les images sont souvent orientés vers autre chose que le Nord. Or, nous le verrons ensuite, il est plus facile d’avoir les images orientées au Nord avant de commencer le géoréférencement.

Nous allons donc réorienter les images au Nord avec Gimp, après avoir préalablement installé le greffon qui lui permet d’écrire sur du Jpeg2000. Le greffon pour Gimp est téléchargable ici

image1924_non_calée

Un fois l’image bien orientée nous allons pouvoir commencer le géoréférencement

 


3- Géoréférencer les images aériennes

Bien !! nous allons pouvoir géoréférencer les images.

Géoréférencer veut dire que nous allons caler les images sur un carte. Le but est d’affecter une référence spatiale, dans une projection géographique donnée, à une image qui n’en a pas. Le résultat obtenu est un Raster, c’est à dire une image divisée en petits rectangles (pixels) qui la décrivent.

image raster : à chaque petit rectangle (pixel) sont associées une ou plusieurs valeurs décrivant les caractéristiques de l’espace

C’est ce raster que l’on peut utiliser dans un  SIG.

Pour cela Qgis a un outil qui s’appelle le géoréférenceur.

Le but est simple mais fastidieux : il s’agit de faire correspondre des points de calages de l’image à caler (des .gcp)  avec une image déjà calée (celle de 2012 dans notre cas).

georeferencer

Le géoréférencement de l’image de 1924 (à droite) avec celles d’aujourd’hui (à gauche) par la définition de points communs de référence (gcp) entre les deux images

 

Une fois les points de calages définis et enregistrés. On choisit le type de transformation de l’image (la plus utilisée est la polynomiale 2).

Qgis va donc plaquer l’image par rotation/étirement, etc. sur la carte actuelle.

 

image_callée

résultat du calage d’un série d’images

 

A savoir :  Qgis génère un rapport de  géoréférencement qui permet de contrôler son travail.

En effet les déformations entre l’angle de prise de vue de l’image de 1924 et une image orthorectifiée peut être important. plus de 10 points par images sont nécessaire pour avoir un résultat convenable

rapport_georeferencement

Rapport de géoréférencement de qgis. on y voit les déformations appliquées à l’image pour la caler

 


4- fusionner les rasters entre eux

Une fois les rasters obtenu il me faut les fusionner ensemble.

Pour avoir le moins de variations d’angles et surtout de colorimétrie, j’ai choisi de les fusionner dans le sens des vols originaux.

 

trajet_vol

trajet des vols en 1924. il donnent le sens de la fusion des raster

L’outil de fusion de Qgis (qui utilise lui même la fonction gdal_merge de la bibliothèque gdal) permet ensuite de fusionner les raster entre eux

gdal_merge

la fonction gdal_merge de Gdal permet de fusionner des raster entre eux

 

 


5- Construction des tuiles

Pour être utilisé sur le web, notre raster doit maintenant être tuilé.

Le tuilage est le re-découpage du raster en petite images à des échelles prédéfinies pour être utilisé par des outils web.

Vous avez sans doute remarqué ça lors de la navigation sur google maps ou sur le géoportail. En effet certaines tuiles apparaissent après les autres ce qui nous permet de voir qu’il ne s’agit pas d’une image unique, trop consommatrice de ressources.

Le tuilage va donc découper le raster sur plusieurs échelles (pyramidage) en autant de petites images (tuilage), dont le nom du dossier et celui de l’images sont les coordonnées x et y dans le système de projection mercator

 

tuilage

Principe du tuilage

pyramidage

principe du pyramidage

Nous allons utiliser le script python gdal2tiles.py qui permet de faire cela en ligne de commande


6 – Développement de l’application avec Leaflet

Une fois les tuiles copiées en ftp sur un serveur on va les appeler avec la librairie javascript Leaflet, utilisé notamment par OpenStreetMap ou Mapbox.

Cette librairie a le mérite d’être assez simple d’emploi et assez puissante. Mais d’autres bibliothèques javascript, comme OpenLayer, existent pour faire de la cartographie web

Pour appeler une carte par exemple ,on insère dans la page html :

1- une<div> html nommée « map »

2- puis  un script javascript avec les balises <script> </script>

on créer ainsi la variable map :

puis ensuite on va intégrer des fonctions que l’on va afficher dans la carte avec la fonction addTo(map)

avec nos tuiles ça donne :

Attention, Gdal2Tiles créer les Tuiles au format TMS, il faut inverser les x et y

On va enfin utiliser un overlay entre les couches qui sera contrôlé par une réglette et généré par la fonction clip();

voici le code entier

 


7- Résultat

Et voilà c’est prêt ! utilisez la réglette en haut de l’image pour passer de 2012 à 1924 !

 

Tagged , , , ,