Category: Géo-traitements

Ressources de traitements de données spatiales

Traiter les dalles SRTM avec GDAL

Cette procédure décrit la méthode pour télécharger et prétraiter les dalles SRTM v4.

Téléchargement :

se connecter sur le serveur CGIAR dans la partie SRTM Data Search and Download et choisir une ou plusieurs dalles de 5°x5°

on obtient une archive contenant 4 fichiers (ex. readme.txt ; srtm_41_05.hdr ; srtm_41_05.tif ; srtm_41_05.tfw)

Mosaïque :

Déplacer tous les fichiers dans un seul répertoire.
Puis lancer la commande :
gdal_merge.py -n -32768 -o out_filename input_files

- n -32768 : pour fixer les valeurs nodata
out_filename : nom du fichier en sorite
input_files : liste de fichiers en entrée

Plus d’options GDAL sur http://www.gdal.org/index.html

Découpage d’une zone

gdal_translate -projwin ulx uly llx lly) src_dataset dst_dataset
ulx uly llx lly sont les coordonnées projetées de la zone à découper, ex. en WGS84 sur la Crête 23 36 27 34
src_dataset : fichier en entrée
dst_dataset : fichier en sortie

Plus d’options GDAL sur http://www.gdal.org/index.html

Ombrage (ou estompage)

gdaldem hillshade -s 111120 input_dem output_hillshade

- s 111120 : ratio pour les unités verticales dans le cas où les coordonnées sont en degrés (e.g. Lat/Long WGS84) et les altitudes en mètres.

Plus d’options GDAL sur http://www.gdal.org/index.html

Dallage avec GDAL

Générer des tuiles avec gdal2tiles.py

création d’un tuilage avec cet utilitaire gdal ici

Description d’une procédure complète pour créer un tuilage et le mettre au format KML de Google Earth ici

Re-tuiler des tuiles existantes

utilitaire ici

exemple pour générer les tuiles en JPEG, compressées à un taux de 80, à partir d’une image ECW

gdal_retile.py -of JPEG -ps 2000 2000 -s_srs EPSG:27572 --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL --config JPEG_QUALITY_OVERVIEW 80 -tileIndex TI_ortho_2003.shp -tileIndexField TileName -targetDir /tmp/tileortho/ /mnt/srvstk0_geomor/Data_SIG/OHM_Mine/Fonds_Reference/Ortho_2003/Ortho_2003_L2e.ecw

voir la documentation sur le site GDAL (rubrique Utilities)

Gestion des dalles raster avec GDAL

Ce document décrit différents scripts basés sur GDAL en Shell pour des dalles raster.
L’exemple d’utilisation est la récupération de données auprès de l’IGN, au format TIFF ou ECW

 

Afficher l’information sur une image raster

gdalinfo src_dataset

Affiche l’information sur un jeu de données raster. Permet notamment d’afficher la référence spatiale et l’étendue des données, le nombre de pixel, la taille du pixel. Documentation officielle gdalinfo

Affecter une référence spatiale à un ensemble de dalles

gdal_translate -a_srs srs_def src_dataset dst_dataset

Affecte le système de coordonnées « srs_def » à l’image en sortie. On utilisera la syntaxe EPSG:n (ex. EPSG :2154 pour du Lambert93 RGF93, 27572 pour du NTF Lambert II Etendu). ex. pour du NTF- Lambert II Etendu :

gdal_translate -a_srs EPSG:27572 im_src.tif im_dst.tif

Script pour un lot de dalles :

#!/bin/bash
for FILE in *.tif
do
BASE=`basename $FILE .tif`
NEWFILE=${BASE}_c.tif
gdal_translate -of GTiff -a_srs EPSG:2154 $FILE $NEWFILE
done

options de gdal_translate

  • compresser une image : l’option -co COMPRESS=DEFLATE permet de compresser l’image TIFF sans perte de qualité. L’option -co TILED=YES permet de spécifier que l’on fait un tuilage.
  • codage de l’image : l’option -co NBITS=1 permet de coder l’image en sortie à 1 bit. Documentation officielle gdal_translate

    Création d’une mosaique au format « Virtual Raster » VRT

    gdalbuildvrt mosaic.vrt *.tif

Crée un mosaïque « mosaic.vrt » à partir d’un ensemble d’une collection d’images tif. Ce format vrt est très léger et semblable aux algorithmes d’Er-Mapper, à savoir que les données d’origine sont conservées et inchangées. Le vrt est simplement un catalogue d’images. Toutes les images de la collection doivent avoir le même système de coordonnées. Pour éviter les zones noires là où il n’y a pas de données on rajoutera une option « -addalpha ». De cette manière, avec un éditeur d’image gérant le canal alpha (ex. ArcMap), on aura de la transparence dans les zones sans données source, et de l’opacité dans les zones avec données source. Ex. :

gdalbuildvrt -addalpha mosaic.vrt *.tif

Documentation officielle gdalbuildvrt

Construire des pyramides

Cette fonction permet d’accélérer l’affichage en créant des aperçus selon les niveaux de zoomes choisis. C’est l’équivalent de la création de pyramides dans ArcGIS.

gdaladdo -ro data_src 2 4 8 16 32

L’option « -ro » permet de créer un aperçu en tant qu’image externe et non inclus dans l’image tiff elle-même. Attention, si vous n’utilisez pas l’option -ro, les pyramides seront construites dans l’image d’origine et elle pourra être corrompue.
Pour compresser le plus possible en JPEG, on utilisera ce type d’options :

gdaladdo -ro --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config JPEG_QUALITY_OVERVIEW 80 data_src 2 4 8 16

Documentation gdaladdo

gdal_merge : mosaiquage

Mosaïque automatique d’un lot d’images. Toutes les images doivent être dans le même système de coordonnées et avoir le même nombre de bandes. Les images peuvent se chevaucher, et avoir des résolutions différentes. Dans les zones de chevauchement, les valeurs prendront celles de la dernière image.

gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*
[-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-v] [-pct]
[-ul_lr ulx uly lrx lry] [-n nodata_value] [-init "value [value...]"]
[-ot datatype] [-createonly] input_files

Exemple pour des dalles SRTM en Geotiff :

gdal_merge.py -o srtm_global.tif *.tif

Aide sur gdal_merge.py

Vecteurs GPS pour ArcGIS

Note préliminaire

ToolBox ArcGis réalisée pour ArcGIS 9.x et 10.0. Non testé avec les versions ultérieures

Construction des vecteurs GPS et des ellipses d’incertitudes.

Boîte à outil ArcGIS

Nécessite en entrée une table contenant au minimum : “Longitude” “Latitude” “Deplacement Est” “Déplacement Nord” ” Incertitude Est” “Incertitude Nord”

INSTALLATION :

- Ajouter la boîte à outils correspondant à votre version dans ArcToolbox

UTILISATION :

- Suivre les étapes 0 à 5 (Si vous disposez déja d’un tableau avec coordonnées projetées, vous pouvez sauter l’étape 1)

- à partir de l’étape 2, utiliser le fichier projeté

Exemple de résultat :

JPEG - 35.5 ko
Vecteurs GPS

TELECHARGEMENT Vecteurs_GPS

BeachBalls pour ArcGIS

Note préliminaire

Macro créée en VB en 2011 avec ArcGis 9.x et 10.0. Non testé avec les versions récentes

Représentation des mécanismes au foyer dans ArcGis

Macro VBA Beach Balls
Macro à installer dans ArcMap

Mécanismes au foyer des séismes en beach balls

Archive à télécharger: Beach_Ball

Programmation VB et VBA dans ArcGIS

Note préliminaire

ceci est quasiment obsolète et je le conserve pour mes propres utilisations.

Un utilisateur qui voudrait scripter ArcGis devra plutôt utiliser Python.

Introduction

L’intégration de VBA dans la version 10 n’est plus assurée. Nous vous conseillons de vous orienter vers des scripts Python pour la Toolbox.

Le document ci-dessous est une introduction à la programmation pour les versions inférieures ou égales à la 9.3.1

Dans de nombreux cas, il est nécessaire de développer soi-même une application qui n’est pas implémentée dans le SIG :

  1. pour une fonction qui n’existe pas en standard ;
  2. pour une opération fastidieuse en standard ;
  3. pour une opération répétitive.

Pour cela, nous utilisons les outils de développement associés à ArcGIS. VBA intégré jusqu’à la version 9.3.1, VB si l’on dispose d’un compilateur. Pour la version 10 et supérieur d’ArcGIS, le support de VBA n’est plus assuré et il faudra préférer utiliser le langage Python, notamment via les scripts ArcToolbox.

Cette documentation est une introduction au développement VBA et VB pour ArcGis.

Doc pdf

documentation et exemples dans ce fichier DOC_CODE_ARCGIS_v0

Filtrage des points sol dans un nuage de points

Différentes approches et outils peuvent être utilisés pour filtrer les points appartenant au sol dans un nuage de points et dont pour créer un MNT

PDAL

PDAL outil Open-source, on pourra voir ce Workshop ou différents exemples sur le site officiel

LASTools

avec LAStools, outil LASGround (payant). C’est la référence en terme de traitements Lidar.

CloudCompare

idées en vrac:

  • utiliser l’algo Canupo (assez difficile à paramétrer).
  • si on dispose d’un MNT Lidar, on peut calculer la distance à cette surface et filtrer grâce à un seuil sur cette distance.
  • On peut aussi utiliser des seuils min et max si la topo est assez plane et qu’on a une idée de la hauteur de la végétation

Agisoft Metashape

La nouvelle version Metashape propose différentes approches pour classifier un nuage de points, voir cet article https://agisoft.freshdesk.com/support/solutions/articles/31000148866-dense-cloud-classification

C’est assez basique comme méthode, mais ça peut aider à faire un pré-traitement.

Divers

post de blog qui référence plusieurs méthodes ici

et méthode et tuto bien écrit avec Grass ici, notamment utilisant l’algorithm MCC

Points2Grid

Points2Grid : A Local Gridding Method for DEM Generation from Lidar Point Cloud Data

site

Points2Grid is used as the default DEM generation algorithm for the OpenTopography point cloud processing system

here is how-to install and use.

Install

sudo apt-get update

sudo apt-get install build-essential gdal-bin libcurl4-gnutls-dev libbz2-dev libboost-all-dev  

  • get the source code : https://github.com/CRREL/points2grid/ . Use either git clone, either download ZIP archive
  • go into the  points2grid directory
  • follow the steps indicated in INSTALL.md, which basically were for me:


sudo mkdir build
cd build
sudo cmake ..        #check no error message in this step and followings
sudo make
sudo make install

  • add to file ~/.bashrc

export PATH=$PATH:/usr/local/bin
LD_LIBRARY_PATH=/usr/local/lib
export LD_LIBRARY_PATH

  • If no error message, you’re done. Test it by running inside the build directory the command:

./points2grid --help

This shoud write down the full help message

Use

  • Use las point cloud as input, if you have LAZ convert it with for example with Lastools:
    LAStools\bin\las2las -i .\20200311_MA_GCP12.laz -o .\20200311_MA_GCP12.las
  • Gridding example command :

./points2grid -i  /path_to_input/20200311_MA_GCP12.las  -o /path_to_output -r 0.0707 --input_format las --all --resolution 0.1 --fill

    • options :
      -r : search radius, which usually is calculated as : resolution*sqrt(2)/2
      –fill : fill empty cells using a 3×3 window
      –all : all the possible calculation are made (see below)

as output, one can get, Oh Miracle!, the following calculation in asc or grid format :

–min                             the Zmin values are stored
–max                             the Zmax values are stored
–mean                            the Zmean values are stored
–idw                             the Zidw values are stored
–std                             the Zstd values are stored
–den                             the density values are stored
–all                             all the values are stored (default)