geotrucs

Gigographie - tartographie

[geotrucs]

pgAdmin 4 - Nouvelle version 3.3

Cette nouvelle version intègre un outil de visualisation des données spatiales.

Quand vous consultez des données ou exécutez une requête, il apparaît maintenant dans le nom de la colonne géométrique un bouton bleu qui permet d’afficher les géométries dans une carte.

Attention, pour avoir un fond de plan OpenStreetMap disponible il faut que vos données soient en WGS84 (4326). Si ce n’est pas le cas faites une reprojection avec la fonction st_transform(geom, srid) de Postgis.

Plusieurs fonds de plan sont disponibles.

[geotrucs]

arcgis : jointures spatiales / listes

Lorsqu’on fait une jointure spatiale sur arcgis, on a deux solutions:

- one to many : renvoie un objet résultat pour chaque match. La jointure peut produire des superpositions de clones géométriques

- one to one : renvoie un seul objet résultat, les données étant agrégées par des fonctions d’agrégation au choix (premier, dernier, somme, moyenne...)

En one to one il existe un moyen d’obtenir une liste des valeurs des objets initiaux en choisissant la fonction d’agrégation joindre et en ajoutant un séparateur. On obtient une chaîne de caractères qu’il est facile de parser en python pour la faire se comporter comme une liste.

Plus de détails ici: https://esriaustraliatechblog.wordpress.com/2015/06/22/spatial-joins-hidden-trick-or-how-to-transfer-attribute-values-in-a-one-to-many-relationship/

[geotrucs]

Compacité d’un polygone - suite

cf ce poste de 2017

On proposait deux indicateurs descriptifs de la compacité d’un polygone:

- son quotient d’isopérimétrie

- sa largeur standard

Comment mettre tout ça en oeuvre en SIG ?

Avant tout, calculer le quotient d’isopérimétrie Qi:

A partir de là, deux cas:

- soit Qi < : dans ce cas le polygone est plus compact qu’un carré (le plus compact des quadrilatères) et la notion de largeur standard n’a pas de sens

- soit Qi > : dans ce cas le polygone est moins compact qu’un carré et dans ce cas sa forme peut être rapportée à une longueur/largeur standards (celles du rectangle de même aire et périmètre que notre polygone)

[geotrucs]

geopackage 101

le problème

Pour ranger les données géographiques d’un projet, on a le choix entre les fichiers shp et une structure plus moderne type base de données. Dans le monde esri coexistent historiquement deux formats: gdb (geodatabase fichier) qui est un format propriétaire esri et mdb (geodatabase personnelle) qui est au départ le format ms access.

Le format gdb est idéal tant qu’on reste dans le monde esri, il reste accessible dans les versions récentes de qgis (via gdal/ogr) moyennant certaines limitations. Par contre, impossible d’accéder aux données attributaires avec access, encore moins excel ou libre office.

Le format mdb est interopérable entre arcgis, qgis et access/libre office mais son gros défaut est que c’est un format obsolète qui n’est plus supporté par microsoft (son inventeur), qu’il comporte des limitations en termes de taille et qu’il sature rapidement la RAM en cas d’utilisation intensive.

Le problème, c’est donc de trouver un format moderne et largement diffusé qui combine les avantages de performance et d’ergonomie du format gdb, et l’interopérabilité du format mdb.

la solution: le geopackage

Le géopackage est un format de données géographique défini par l’OGC (donc libre et ouvert) basé sur le format de base de données sqlite.

Un fichier sqlite contient à la fois une base de données et le moteur de base de données qui permet d’exploiter la donnée, le tout totalement compatible SQL. Le format sqlite est très discret mais universellement répandu grâce à sa capacité à fonctionner sous n’importe quelle infrastructure matérielle et n’importe quel OS. On le trouve dans firefox, chrome, android... C’est donc le format interopérable par excellence.

Un geopackage c’est simplement une base sqlite qui contient une cartouche spatiale permettant de gérer des données géographiques vectorielles. La donnée attributaire attachée reste simplement au format sqlite classique (même principe que la donnée attributaire dans une mdb).

Pour plus de détails:

http://www.geopackage.org/

https://www.sqlite.org/index.html

https://en.wikipedia.org/wiki/GeoPackage

https://fr.wikipedia.org/wiki/SQLite

Dans la suite de cette note, on démontrera toute la chaîne d’opérations depuis la création jusqu’à l’édition de données dans un geopackage avec tous les outils à notre disposition: arcgis pro/desktop, qgis 3, access et libre office.

Prérequis

arcgis desktop/pro: aucun prérequis, les outils sont disponibles en standard

qgis: idem en version 3, interface presse-bouton

ms access: nécessite le driver sqlite 3.0 pour accéder à la donnée attributaire (http://www.ch-werner.de/sqliteodbc/)

libre office: idem, installation du même driver

Créer un geopackage

avec arcgis desktop/pro: data management tools/workspace/create SQLite Database. Bien spécifier le type spatial: GEOPACKAGE_1.2

avec qgis 3: couche>gestionnaire de sources de données>geopackage>nouveau (et plusieurs autres endroits dans l’interface)

avec access/libre office, on créera une base sqlite, mais sans cartouche spatiale

lire/écrire des données géo dans un geopackage

avec arcgis desktop/pro: le geopackage est un workspace comme n’importe quel autre. On peut enregistrer une couche, créer une nouvelle couche, utiliser le geoprocessing copy features exactement comme avec une gdb. Attention: l’écriture dans un gpkg nécessite une licence basic

avec qgis 3, c’est le même principe. Le geopackage est une source comme une autre, on la manipule exactement de la même manière

avec access: on accède aux données attributaires en se connectant au geopackage comme à une source externe. Attention toutefois: les tables avec un prefixe gpkg constituent la cartouche spatiale. Il ne faut surtout pas y toucher, sous peine de casser le fichier.

[geotrucs]

Altitude en open data

De plus en plus de données d’altitude sont disponibles en accès libre. Les trois sources citées ci-dessous ont en commun leur couverture mondiale, leur gratuité et leur précision. Pour autant ces données ne sont pas complètement similaires et bien les connaître permet de s’orienter pour faire un choix.

Réponse courte: JAXA ALOS !

Réponse longue: voir ci-dessous

Space Shuttle Radar Topography Mission (SRTM)

Produit par la Nasa en 2000, la résolution est constante: 30 m sur toute la couverture (mondiale) depuis 2015. Le capteur est un radar inSAR (https://en.wikipedia.org/wiki/Interferometric_synthetic-aperture_radar). Ce qu’il faut savoir de ce capteur, c’est qu’il peine dans les zones de fort relief en bordure de scène: effets de distorsions ou zones complètement noires.

Téléchargement en dalles de 1 deg ici: http://earthexplorer.usgs.gov/ (sur inscription gratuite)

ASTER Global Digital Elevation Model

NASA - Japon: couverture mondiale avec résolution de base de 30 m. La technologie utilisée est une corrélation stéréoscopique d’images satellitaires dans le visible et le proche infrarouge. Une V2 sortie en 2011 corrige beaucoup de défauts de jeunesse bien que les images présentent toujours beaucoup de bruit. Ce bruit est gênant en zone de plaine, mais ASTER est considéré comme plus fiable que SRTM en zone accidentée.

Téléchargement en dalles de 1 deg ici: http://reverb.echo.nasa.gov/reverb/

JAXA’s Global ALOS 3D World

MNT mondial à 30 m dont la dernière version est sortie en mars 2017 comme une version dégradée de AW3D, un MNT mondial à 5 m (payant bien sûr). ALUS n’apporte pas d’amélioration de résolution par rapport à ses deux prédécesseurs, mais le gain de qualité est incontestable.

La donnée est disponible en dalles de 5 deg ici : http://www.eorc.jaxa.jp/ALOS/en/aw3d30/ (inscription gratuite)

[geotrucs]

arcgis pro 2.0 (2/5) le mode projet

Fini les mxd !

Et c’est sûrement le changement le plus important qu’implique l’utilisation d’arcgis pro en production.

Au démarrage d’arcgis pro, on choisit un projet (nouveau ou en cours) et toute la session se déroulera dans ce projet. La doc officielle sur la notion de projet dans pro: https://pro.arcgis.com/fr/pro-app/help/projects/what-is-a-project.htm

En bref: qu’est-ce qui change par rapport aux mxd

Un mxd, comme un wor (mapinfo), contenait tous les éléments nécessaires à la production d’une carte papier:

une mise en page avec ses éléments graphiques

un, ou plusieurs blocs carte, qui s’insèrent dans cette mise en page

des appels à un ensemble de données (shp, jointures...) ainsi que des définitions de symbologies

Point de détail, mais important, un mxd était un fichier binaire qui n’était modifiable (en dehors d’arcgis) qu’au travers d’une api arcpy. Pour faire simple, un mxd = une carte papier. Dans une étude, on multiplie les mxd.

Un aprx (projet arcgis pro) c’est bien plus que ça et c’est conçu pour coller à la notion de projet sig avec des cartes papier, des cartes web, des métadonnées, des process, des tâches... Les éléments du projet sont regroupés dans le panneau “catalogue”:

mises en page (au pluriel)

cartes

boîtes à outils

bases de données (locales et connexions)

dossiers

styles

géocodeurs

tâches

De nouveaux éléments font leur apparition: des tâches (guides détaillés de process), des notifications, des boîtes à outils. L’objectif est de favoriser la collaboration, la mutualisation.

Techniquement, comment ça se présente

Un dossier projet (ici __nada__ qui est mon projet par défaut) contient une gdb, une boîte à outils, un fichier aprx (définition de projet), un dossier de log d’imports et un historique des géotraitements (.pyHistory).

Lorsqu’on ajoute une donnée distante dans le projet, on ajoute au projet une connexion à une gdb, à un dossier et ces connexions sont référencées. On peut ajouter des dossiers à l’arborescence de base du projet (ex: pdf, jpg, sources...).

On voit qu’on s’achemine vers une normalisation des dossiers étude (un de mes dadas !).

Un détail, mais non des moindres : le fichier aprx est en fait un simple jeu de fichiers xml zippés qui contient toute les métadonnées et les paramètres du projet. C’est complexe, mais il reste possible de le lire/modifier avec un simple éditeur de texte.

[geotrucs]

arcgis pro 2.0 (1/5) frais style

Une interface sombre, comme nos collègues infographes, anyone ?

Les outils à gauche, la carte plein écran à droite ? Avec un peu de place pour outlook ou un google keep ?

Un bandeau comme office 2013 et des panneaux ajustables, superposables à volonté... De l’apparence, bien futile tout ça, mais pas désagréable, d’autant plus que l’interface est enfin fluide: pas comme arcgis desktop !

[geotrucs]

Géométries doublons avec arcpy

...et pandas

Le principe est d’effectuer un “group by” sur l’objet géométrique et d’appliquer une fonction d’aggrégation: pour celà on peut passer par le token SHAPE@WKT de la fonction arcpy.da.SearchCursor, puis de réaliser une agrégation avec pandas.

Un exemple:

Ici le but est de ne garder que l'objet dont la valeur 'struct_num' est la plus faible

En fin de parcours on extrait les OBJECTID correspondants et on fabrique un querydef

[geotrucs]

Installer proprement qgis sur windows

Où télécharger les installateurs?

Ici: http://www.qgis.org/fr/site/forusers/download.html

32 ou 64 bits?

64 bits !

En 2017 tous les processeurs sont multicoeurs. On choisira donc par défaut un installateur 64 bits, sauf contexte très particulier (machine virtuelle dotée d'un seul coeur et d'un OS 32 bits par exemple)

installateur indépendant ou réseau?

Installateur réseau !

Pour une installation simple sur un poste de travail on a spontanément tendance à choisir l'installateur indépendant. A mon avis, dans un cadre professionnel, il vaut mieux choisir l'installateur réseau qui fonctionne aussi très bien sur un poste indépendant et qui offre en plus des possibilités d'installation de dépendances, de logiciels complémentaires et surtout de mises à jour très intéressantes.

LTR ou latest ?

Lastest !

La version LTR actuelle est la 2.14, tandis que la dernière version stable est la 2.18.9. La LTR est réputée particulièrement stable et conseillée dans le cadre d'une utilisation professionnelle en production. Cependant le développement de qgis est rapide et de nombreuses nouvelles fonctionnalités sont régulièrement ajoutées dans la version à jour. Pour moi le bénéfice - risque penche largement en faveur de la dernière version, d'autant plus que les bugs rédhibitoires restent rares.

Comment installer qgis avec l'installateur réseau ?

télécharger l'installateur et le lancer, de préférence en tant qu'admin

choisir express desktop install

laisser tout par défaut et accepter les licences proposées

L'installateur installe qgis, grass avec son plugin, saga-gis et un certain nombre de dépendances, notamment gdal-fileGDB qui offre un acces rw aux géodatabases fichier d'esri.

Comment mettre à jour qgis avec l'installateur réseau ?

relancer l'installateur (dans le dossier SSGeo4W) et choisir l'installation avancée, toujours en tant qu'admin

accepter les nouvelles licences et tout laisser par défaut. On peut aussi rajouter des dépendances qui n'auraient pas encore été installées

[geotrucs]

Compacité d’un polygone

En analyse spatiale, on est souvent amené à fixer un seuil en dessous duquel une parcelle n’est plus accessible à un usage ou un phénomène (constructibilité, niche écologique...). Le premier critère qui vient à l’esprit, c’est la surface: par exemple

en dessous de 800 m² une parcelle n’est pas constructible

Mais la superficie n'est pas le seul critère: une bande continue de 100 km par 3 m le long d'une autoroute ne fait pas un habitat utile pour un lynx.

Vient donc un deuxième critère qu'on pourrait définir comme la compacité: intuitivement un rapport périmètre/aire. Malheureusement celui-ci varie en fonction de la taille de l'objet: on préfèrera

q = 4piA/P²

qui est indépendant de la taille de l'objet:

pour un disque q = 1

pour un triangle équilatéral q =

pour un carré q = 1/4pi

pour un rectange q = 4piA/P²

q est en fait le quotient isopérimétrique: https://fr.wikipedia.org/wiki/Isop%C3%A9rim%C3%A9trie

Avec un couple de critères A et q, on peut éliminer facilement les polygones les moins compacts, soit ceux dont le rapport du périmètre à l'aire sont les plus éloignés du disque.

Mais ce quotient n'est pas très parlant. On peut en dériver un indicateur plus beaucoup plus explicite: la "largeur standard".

Assimilons le polygone à un rectangle et définissons le couple longueur (L) et largeur standards (W) tels que L et W vérifient

P = 2(L + W)

A = LW

La résolution du polynôme du second degré passe par la recherche du discriminant (Attention programme de math en seconde) :

Δ = 1/4P² - 4A

Rq: si Δ <= 0, le polygone n'est pas assimilable à un rectangle: il n'y a pas de couple L et W qui vérifient les propriétés d'aire et de périmètre énoncées plus haut. Dans ce cas

on en déduit

L = 1/4P - 1/2√Δ

et bien sûr

W = 1/2P-L

On des longueurs, beaucoup plus faciles à comprendre. Par exemple, on peut définir que la largeur standard limite est de 20 m pour une parcelle constructible, en plus du seuil de 800 m², ou que la largeur standard limite est de 5 km pour le territoire d'un lynx.

Le couple aire/largeur standard dérivent directement du quotient isopérimétrique et décrivent de la même manière la compacité du polygone. Cependant ce couple de valeurs est plus facile à appréhender dans le cadre de la définition d'un seuil.

[geotrucs]

in_memory

C’est là qu’intervient l’espace de stockage temporaire in_memory: plutôt que d’enregistrer la donnée en sortie sur le disque, on la stocke dans la RAM.

avantages: c'est temporaire! de plus les I/O sont infiniment plus rapides

inconvénients: c'est temporaire! si on a chaîné plusieurs opérations et que le traitement plante en cours de route, il y a de fortes chances que les données intermédiaires soient perdues

Comment utiliser in_memory: dans les boîtes de dialogue de géotraitement, il suffit de remplacer les chemins

X:\chemin\vers\la\base.gdb\donnee_en_sortie

par

in_memory\donnee_en_sortie

Idem dans appel à une fonction arcpy.

La donnée produite est très similaire à une donnée gdb à quelques détails près (pas de champs shape_area, shape_length...)

Plus d'infos:

http://desktop.arcgis.com/fr/arcmap/latest/analyze/modelbuilder/the-in-memory-workspace.htm

[geotrucs]

A propos du multithread/multiprocess

Pourquoi 25% ? Tout simplement parce que nos CPU “modernes” (enfin depuis bientôt 10 ans) ont au minimum 4 coeurs et parce que les algos que nous utilisons n’en tirent pas parti. Lorsque nous lançons un traitement, celui-ci est affecté à un “process”, lui-même affecté à un coeur. Les opérations sont donc mises en file et effectuées l’une après l’autre. Pendant ce temps, un autre coeur est affecté à faire tourner le système (OS) pour moins de 10% de charge et les deux autres se tournent les pouces. Enorme perte d’efficacité!

La seule manière d’optimiser ça, c’est de répartir le travail sur plusieurs coeurs: dans le cas d’un quad-core, si on peut entièrement dédier la machine au traitement, on peut mettre les 4 coeurs en charge, mais ça sera au prix d’un gros ralentissement du système: inutile d’espérer écrire un mail pendant le traitement. Dans ce cas on réservera 3 coeurs au traitement et 1 à la bureautique.

Comment faire?

L’approche générale est toujours en 3 temps:

scinder la donnée en autant de jeux qu’on veut faire travailler de coeurs

traiter ces jeux simultanément par autant d’instances de l’algorithme de de traitement

rassembler le résultat

C'est très efficace dans un traitement qui peut être décomposé ligne par ligne (la plupart des géotraitements), impossible à mettre en oeuvre pour une requêtre SQL. Dans ce cas-là mieux vaut compter sur le moteur SGBD du serveur qui est optimisé pour aller vite.

Pour faire du multiprocessing, on peut utiliser des bibliothèques dédiées (multiprocessing pour python) qui reproduisent exactement ce schéma, mais l’approche “à la main”, qui ne nécessite aucune connaissance en programmation, est souvent très efficace... Avec arcgis c'est le moment de tirer de l'espace de stockage r'in_memory' pour stocker les données intermédiaires !

[geotrucs]

multiprocessing sur grass

facile!

A faire tourner dans une session grass:

[geotrucs]

Ma source est-elle une gdb ?

Pourquoi se poser la question ?

...check doc...

Comment vérifier ?

[geotrucs]

Classe d’entités vers DataFrame pandas

Rapide et crado : on passe par un dictionnaire qui ralentit beaucoup la manip.

[geotrucs]

Clustering : techniques de typologie

http://www.sthda.com/english/wiki/hcpc-hierarchical-clustering-on-principal-components-hybrid-approach-2-2-unsupervised-machine-learning#compute-hierarchical-clustering-on-the-pca-results

Une approche extrêmement synthétique des techniques de clustering avec R et sa librairie d’analyse multivariée FactoMineR. Les techniques de clustering sont particulièrement utiles pour effectuer des typologies sur facteurs multiples.

Deux cas sont explorés:

variables continues

variables discrètes

Les résultats sont affichés à l’aide de la librairie factoextra, spécialisée dans les graphes de clusters.

[geotrucs]

Cascade maps

Comment décomposer une carte complexe en un récit de cartes en cascade?

Une application par l’exemple avec l’épisode tragique de l’invasion des hamsters zombies: https://nation.maps.arcgis.com/apps/Cascade/index.html?appid=954145df6cf84e2d8bbea996438c99fb

Et un autre exemple de présentation d’une étude complète sans la moindre carte...

https://canalrivertrust.maps.arcgis.com/apps/Cascade/index.html?appid=f9347d20313f48e2a0d689de3f4998d3