https://grasswiki.osgeo.org/wiki/Hydrological_Sciences
https://grasswiki.osgeo.org/wiki/Creating_watersheds
https://baharmon.github.io/hydrology-in-grass
SIGéo CEREGE
Follow this tutorial https://www.youtube.com/watch?v=-xzoVF7Z7u0
Here is a Qgis style layer as example :
Differences between ellipsoidal height and orthometric altitude (or geoid height) are explained here https://spatialthoughts.com/2019/10/26/convert-between-orthometric-and-ellipsoidal-elevations-using-gdal/
The definition of the projection is the PROJ4 line that can be found on https://spatialreference.org/
gdalwarp -s_srs "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs" -t_srs "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +geoidgrids=egm96_15.gtx" input.tif output.tif
gdalwarp -s_srs "+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs +type=crs +geoidgrids=RAF20.gtx" -t_srs "+proj=longlat +datum=WGS84 +no_defs +type=crs" input.tif output.tif
gdalwarp -s_srs "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs
+geoidgrids=egm96_15.gtx"
-t_srs "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs" input.tif output.tif
you can choose your input and output spatial reference using the proj4 syntax https://epsg.io/32633 , and note the add of the geoid file with +geoidgrids=egm96_15.gtx in the command.
if you want to use EGM2008 geoid, get the file from https://github.com/OSGeo/proj-datumgrid/blob/master/world/egm08_25.gtx and copy it to the proj directory. In Linux it’s
/usr/share/proj
then you can simply change the geoid file in the command +geoidgrids=egm08_25.gtx
Lors de l’écriture d’un raster, surtout en THR (très haute résolution), il est important de compresser au maximum les images pour avoir des fichiers qui soient gérables. Il est important d’avoir en tête que certaines compression se font avec une perte de qualité, et d’autres sans perte. Si l’on veut simplement faire de la visualisation, on peut se permettre un peu de perte, par contre pour du traitement ou pour des données type MNT, on ne voudra aucune perte d’information.
Toutes ces histoires de compression sont bien expliquées sur ce blog http://blog.cleverelephant.ca/2015/02/geotiff-compression-for-dummies.html
Possibilité de tuiler (daller) pour permettre de gérer des trop grosses images. C’est parfois indispensable, par exemple si l’on passe au dela de 20 à 30000 pixels de large pour une image.
On pourra consulter pour plus d’infos cette page https://sigeo.cerege.fr/?p=116
le GeoTIFF permet la gestion de très grosses images avec un système de tuilage interne. Tout est donc stocké dans 1 fichier TIFF.
Pour les images >= 4Go, il faudra les écrire en BigTiff.
Parfois cela est réalisé automatiquement par l’application qui va écrire l’image. Si l’on utilise GDAL, il est préférable de le signifier explicitement avec l’option “-co BIGTIFF=YES”
pour forcer l’écriture des tuiles en interne, c’est avec l’option -co TILED=YES
exemple de commande avec ces 2 options: gdal_translate im_in.tif im_out.tif -co BIGTIFF=YES -co TILED=YES
par exemple pour un MNT ou de l’imagerie que l’on veut garder intacte, j’utilise la compression avec DEFLATE
ex . gdal_translate dem_in.tif dem_out.tif -co COMPRESS=DEFLATE
on peut utiliser la compression JPEG pour réduire efficacement la taille des images multi-canaux en 8 bits (ne fonctionne pas pour les images 16 ou 32 bits) :
gdal_translate \
-co COMPRESS=JPEG \
-co PHOTOMETRIC=YCBCR \
-co TILED=YES \
input.tif output.tif
Le stockage en JPEG est encore plus performant dans l’espace de couleurs YCBCR.
Pour des images en 16 ou 32 bits, on utilisera LZW ou PACKBITS pour la compression avec perte de qualité.
Les pyramides permettent d’accélérer l’affichage en créant des aperçus selon les niveaux de zoomes choisis. Ici la perte de qualité n’est pas un problème car c’est simplement pour l’affichage, et la donnée originale n’est pas altérée.
Les pyramides peuvent être stockées dans un fichier séparé (format OVR), ce qui est en général préférable pour plus de clarté, ou en interne dans un GeoTiff.
exemple de pyramides pour de l’imagerie 8 bits :
gdaladdo \
--config COMPRESS_OVERVIEW JPEG \
--config PHOTOMETRIC_OVERVIEW YCBCR \
--config INTERLEAVE_OVERVIEW PIXEL \
-r average \
-ro \
im_input.tif \
2 4 8 16 32
pour un MNT (16 ou 32 bits en général), on fera plutôt avec du DEFLATE ou du LZW
gdaladdo -ro --config COMPRESS_OVERVIEW LZW dem_input.tif 2 4 8 16 32 64
Compression pour dummies : http://blog.cleverelephant.ca/2015/02/geotiff-compression-for-dummies.html
Format GeoTifff : https://gdal.gloobe.org/gdal/formats/gtiff.html
Format Geotiff : https://www.gdal.org/frmt_gtiff.html
gdaladdo pyramides : https://www.gdal.org/gdaladdo.html
remarque générale: les fichiers rasters intermédiaires produits peuvent être écrits en format VRT (raster virtuel). Cela permet d’avoir de tous petits fichiers.
exemple avec un fichier in.tif en entrée, et écriture d’un fichier out.vrt en sortie
gdal_translate in.tif out.vrt -ot byte -of VRT -co NBITS=1
l’option -co NBITS=1 permet de passer l’image en binaire
pour écrire un fichier en sortie au format VRT, utiliser l’option -of VRT
gdalwarp -dstalpha out.vrt out1.vrt -co NBITS=1
Le canal alpha est créé avec toutes les données en nodata grâce à l’option -dstalpha
on peut aussi prendre en compte une valeur représentant le nodata dans le fichier en entrée grâce à l’option -srcnodata value
Bien vérifier le raster out1.vrt en sortie pour voir dans quel canal se trouve la bande alpha à utiliser par la suite.
gdal_polygonize.py out1.vrt -b 1 -f "ESRI Shapefile" out2.shp
l’option -b 1 permet de prendre ce qui est dans la bande 1
on calcule un raster avec une expression logique là où il y a des valeurs
gdal_calc.py -A input --outfile=tmp.tif --type=Byte --calc="isnan(A)"
on convertit en polygone le raster:
gdal_polygonize.py tmp.tif -f "ESRI Shapefile" output
On suppose dans cette procédure que nos données constituent un maillage régulier et qu’il n’est pas nécessaire de faire d’interpolation. Le fichier ASCII XYZ est composé de trois colonnes X, Y et Z, ex.
880000 | 6340000 | 828 |
880025 | 6340000 | 831 |
880050 | 6340000 | 833 |
880075 | 6340000 | 835 |
1-Changer l’extension du fichier en .txt ou .csv
2- Modifier le fichier schema.ini
Plus d’infos sur ce fichier schema.ini à http://www.tectonics.caltech.edu/gi…
Ex. de fichier schema.ini correspondant au fichier input.csv ci-dessus : [input.csv]
ColNameHeader=False
Format=Delimited( )
Col1=X Double
Col2=Y Double
Col3=Z Double
Si le fichier n’est pas bien reconnu par ArcCatalog, il est possible que ce soit un problème de codage des caractères. On pourra le modifier en ajoutant la ligne suivante dans schema.ini :
CharacterSet=UNICODE
Changer UNICODE par le codage correct (ANSI, autre …)
3- Vérifier la reconnaissance du fichier par ArcCatalog
4- Créer le shapefile de points
5-Convertir le shapefile de points en raster
ex. conversion d’un fichier XYZ ASCII en TIF et assignation de la référence spatiale RGF93-Lambert93 (Code EPSG 2154)
gdal_translate -a_srs EPSG:2154 input.XYZ output.tif
en batch on pourra faire un script du type :
#!/bin/bash
#Ce script assigne le code EPSG saisi en paramètre aux images en sortie
# Jules Fleury, SIGeo/CEREGE, mars 2011
# saisir le code epsg en parametre
# Assign_SRS_to_TIF.sh
echo “Input spatial reference EPSG (exemple: for WGS84 enter 4326) : ”
read c_EPSG
for FILE in *.XYZ # modifier ici en cas d'autre format
do
BASE=`basename $FILE .XYZ` # modifier ici en cas d'autre format
NEWFILE=${BASE}_c.tif
gdal_translate -a_srs EPSG:$c_EPSG $FILE $NEWFILE
done
solution prise ici : http://osgeo-org.1560.x6.nabble.com/gdal-dev-Opening-gridded-xyz-data-that-is-out-of-order-td5334341.html
if we obtain this kink of error:
> ERROR 1: Ungridded dataset: At line 4125001, change of Y direction
> gdalinfo failed – unable to open ‘DTM_swissALTI3D_XYZ.txt’.
You must sort the data beforehand
sort -k2 -n -k1 input.xyz -o output_sort.xyz
gdal_translate input.xyz output.tif -a_nodata -9999
On utilise la fonction d’interpolation de GDAL gdal_grid pour ne prendre pour chaque pixel qu’un seul point en entrée. La méthode consiste à premièrement créer un CSV, puis un VRT et utiliser ce fichier comme source de gdal_grid.
ex. gdal_grid -a nearest:radius1=1:radius2=1:nodata=-9999.0 -ot Float32 -a_srs EPSG:2154 -txe 1074000 1079000 -tye 6325000 6330000 -outsize 1000 1000 -l grid grid.vrt grid_out.tif
traitement en batch
1ère étape : conversion des XYZ en CSV avec création de VRT associé
ex. de script
#!/bin/bash
#Lecture de la référence spatiale
echo “Input spatial reference EPSG (exemple: for WGS84 enter 4326, for Lambert93 2154, for L2E 27572) : ”
read c_EPSG
#traitement des fichier en batch
for IFile in *.xyz # modifier ici en cas d’autre format
do
BASE=`basename $IFile .xyz`
OFile=${BASE}.csv
#ajout d’une ligne d’entête “X Y Z”
sed ‘1i\X Y Z’ $IFile > $OFile
#remplacement du caractère espace par une tabulation
TFile=”Temp”
tr ‘ ‘ ‘\t’ < $OFile > $TFile
mv $TFile $OFile
#création du fichier VRT
OFileVRT=${BASE}.vrt
echo "<OGRVRTDataSource>" > $OFileVRT
echo " <OGRVRTLayer name=\"$BASE\">" >> $OFileVRT
echo " <SrcDataSource>$OFile</SrcDataSource>" >> $OFileVRT
echo " <GeometryType>wkbPoint</GeometryType>" >> $OFileVRT
echo " <LayerSRS>EPSG:$c_EPSG</LayerSRS>" >> $OFileVRT
echo " <GeometryField encoding=\"PointFromColumns\" x=\"X\" y=\"Y\" z=\"Z\"/>" >> $OFileVRT
echo " </OGRVRTLayer>" >> $OFileVRT
echo "</OGRVRTDataSource>" >> $OFileVRT
done
2ème étape : conversion des VRT (avec leur CSV) en TIF
on utilise gdal_grid,ex. de script
#!/bin/bash
#Programme pour interpoler des fichiers de points XYZ stockés en CSV/VRT
#selon différents algorithmes et paramètres.
#le fichier grid en sortie est stocké en TIF
#puis rééchantilloné selon la résolution voulue.
#Jules Fleury – SIGéo/CEREGE
#avril 2013
#Lecture de la référence spatiale
echo “Input spatial reference EPSG (exemple: for WGS84 enter 4326, for Lambert93 2154, for L2E 27572) :”
read c_EPSG
#Traitement des fichiers en batch
for IFile in *.vrt # modifier ici en cas d’autre format
do
echo “Processing $IFile”
BASE=`basename $IFile .vrt`
OFile=${BASE}.tif
val_nodata=-9999 #valeur pour nodata
#Traitement du fichier en utilisant l’algorithme Nearest
gdal_grid -a nearest:radius1=0.5:radius2=0.5:nodata=$val_nodata -ot Float32 -a_srs EPSG:$c_EPSG -l $BASE $IFile $OFile
#ré-échantillonage du fichier en sortie
ressize=10 #valeur de ré-échantillonage
OFileRes=${BASE}_$ressize.tif #fichier ré-échantilloné en sortie
gdalwarp -tr $ressize $ressize -r cubic $OFile $OFileRes
echo "Done"
done
On fera attention au séparateur de champs du fichier (virgule par défaut, tabulation, espace) et à la présence ou non d’une entête.
2 possibilités, en utilisant ArcGIS ou en utilisant OGR
1- Renommer votre fichier avec l’extension .csv ou .txt
2- Modifier le fichier schema.ini
Plus d’infos sur ce fichier schema.ini à http://www.tectonics.caltech.edu/gi…
Ex. de fichier schema.ini correspondant au fichier input.csv ci-dessus :
[input.csv]
ColNameHeader=False
Format=Delimited( )
Col1=X Double
Col2=Y Double
Col3=Z Double
Si le fichier n’est pas bien reconnu par ArcCatalog, on peut définir le codage des caractères en ajoutant la ligne suivante dans schema.ini : CharacterSet=UNICODE Changer UNICODE par le codage correct (ANSI, autre …)
3- Vérifier la reconnaissance du fichier par ArcCatalog
Vérifier dans ArcCatalog que votre fichier est bien reconnu. Pour cela votre fichier doit s’afficher dans l’onglet Aperçu. Vérifier aussi que le type des champs est correct. Pour cela faire un clic droit sur le fichier et vérifier l’onglet Champs
4- Créer le shapefile de points
Clic droit sur le fichier – Créer une classe d’entités – A partir d’une table XY – Renseigner les champs contenant les coordonnées X, Y et éventuellement Z ; renseigner aussi le système de coordonnées ; spécifier le fichier en sortie – Cliquer sur OK Vérifier votre shapefile.
Plus d’infos sur http://www.gdal.org/ogr/drv_csv.html et http://gdal.gloobe.org/ogr/formats/…
1- Renommer votre fichier avec l’extension .csv
2- Ajouter la ligne d’entête La ligne d’entête doit contenir les noms des colonnes séparés par le délimiteur. Ex.
X Y Z
869324.598 1783031.225 125.000
869781.445 1783057.845 126.000
869340.038 1782760.764 127.000
Pour cela, on peut soit éditer le fichier avec un éditeur de texte soit utiliser une commande sed dans le shell : Ex. (ajout d’une ligne ’’X Y Z’’ en début de fichier)
sed '1i\X Y Z' input.csv > output.csv
3- Remplacement du caractère Espace par le caractère Tabulation
OGR ne reconnaissant par le caractère Espace comme séparateur, il faut le remplacer par une tabulation ou une virgule. Pour cela on utiliser une commande shell :
tr ' ' '\t' <input.csv >>output.csv
on peut aussi utiliser sed, exemple :
sed 's/[ ]/\t/g' test3.csv > test3_b.csv
4- Lecture d’un CSV contenant des données spatiales
Il est possible d’extraire de l’information spatiale (points) d’un CSV ayant des colonnes X et Y en utilisant le driver VRT (http://www.gdal.org/ogr/drv_vrt.html ) Il faut écrire le fichier VRT associé (input.vrt) au fichier CSV (input.csv)
<OGRVRTDataSource>
<OGRVRTLayer name="input">
<SrcDataSource>input.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>EPSG:27572</LayerSRS>
<GeometryField encoding="PointFromColumns" x="X" y="Y"/>
</OGRVRTLayer>
</OGRVRTDataSource>
et ogrinfo input.vrt renverra :
INFO: Open of `input.vrt'
using driver `VRT' successful.
1: input (Point)
5- Conversion du CSV en Shapefile
ogr2ogr output.shp input.vrt
On peut ajouter un filtre spatial sur les coordonnées en utilisant “-spat”
ogr2ogr -spat 822500 6245000 836500 6264000 -a_srs epsg:2154 Veran_T_total.shp T_total.vrt
Feuilles de calcul (Excel et Calc) pour convertir des coordonnées entre :
Degrés Minutes Secondes ;
Degrés décimaux ;
Degrés Minutes décimales.
Téléchargement DDtoDMS_v2
Cette documentation traite des procédures de traitement de profils de chenal fluvial, partant des profils en travers et jusqu’aux profils en long. Ce sont plusieurs macros VBA sous Excel.
La documentation décrit les macros et inclus le code source : docu_codeexcel
A partir d’un fichier texte contenant les coordonnées de points
Format de fichier en entrée : NPT,nb_de_pts
X,Y,Z
11225,2254,32
12143,2341,11
…
=>Calcule la distance entre chacune des paires de points Format de fichier en sortie
FromPoint ToPoint Dist
0 1 48.8756299192227
0 2 90.1678706635675
0 3 56.1265436313451 …
Archive : Distance_paires