Using Proj transformation

PROJ is a generic coordinate transformation software that transforms geospatial coordinates from one coordinate reference system (CRS) to another. It is currently used in many GIS softwares, but can also be used through an API or by command-line tools.

Below some usefull command lines:

Simple conversion from a known coordinate system (CS) to another known CS :

here from UTM33N-WGS84 to WGS84

echo 357000 4676000 | cs2cs +init=epsg:32633 +to +init=epsg:4326

Should return

13d16'2.189"E 42d13'23.141"N 0.000

More documentation https://proj.org/usage/quickstart.html

to add some formatting to output, specify the format with “-f” as a printf format string. For example “-f %.8f” will return the coordinate in decimal degrees with 8 decimals

13.26727481 42.22309481 0.00000000

to input several coordinates at a time, either use a file and pipe with cat

cat coords_utm.txt | cs2cs +init=epsg:32633 +to +init=epsg:4326 -f %.8f Should return the list of converted coordinates
13.26727481 42.22309481 0.00000000
13.54977927 42.05589846 0.00000000

and to output into a file, redirect with > symbol

cat coords_utm.txt | cs2cs +init=epsg:32633 +to +init=epsg:4326 -f %.8f > coords_dd_wgs84.txt

to input manually several coordinates, use EOF or whatever specific characters you like :

cs2cs EPSG:4326 EPSG:32631 <<EOF
45N 2E
EOF

Coordinate system transformation : ITRF, ETRF, WGS84

Clear introduction about Global and local referential (Lantmateriet) : “Positions determined by the GNSS method Precise Point Positioning (PPP) are in the same reference frame as the orbits, i.e. usually a realization of ITRS, e.g. ITRFyy, IGSyy or WGS84, where “yy” represents the year of the realization. The coordinates change with time in the ITRS realizations, because of the plate tectonics. Hence,the determined coordinates are given in the epoch of the observations. For practical applications like mapping and referencing spatial data , a static system/frame, which does not change with time, is desired. For this purpose, ETRS89 has been developed for Europe. ETRS89 coincides with ITRS at epoch 1989.0.”

Another simple reading about dealing with ITRS, ETRS and WGS84 is at Confluence website.

Nota : the transformations made by GIS software from WGS84 to local referential are precise at 1 meter only. For centimeter accuracy, use a geodetic software taking into account the velocities of the ITRF referential relative to the local referential (e.g. ETRF for Europe).

There are many WGS84 realizations. The latest compares with ITRF08 and ITRF14, see table below:

YearRealization (Epoch)For all practical purposes equivalent to:
1987WGS 1984 (ORIG)NAD83 (1986)
1994WGS84 (G730)ITRF91/92
1997WGS84 (G873)ITRF94/96
2002WGS84 (G1150)ITRF00
2012WGS (G1674)ITRF08
2013WGS (G1762)Compares to ITFR08 within 1cm Root Mean Square (RMS) overall
Table of WGS84 compared to ITRF from https://www.e-education.psu.edu/geog862/node/1804

For US, WGS84 (G1762) is equivalent at 1cm to ITRF14. For France also ITRF14 ~ ITRF08 at less than 1cm, the transformation is given here http://itrf.ensg.ign.fr/trans_para.php

For France, the official referential is RGF93, which in its latest version is defined as ETRF2000 (epoch 2009.0) ( https://geodesie.ign.fr/index.php?page=rgf93 ).

Position and velocities of the IGS stations in ITRF2014 at epoch 2010 is given at this adress http://itrf.ensg.eu/ITRF_solutions/2014/doc/ITRF2014_GNSS.SSC.txt

Tools to Convert between referentials

Reference online tool to convert between ITRF and ETRF : http://www.epncb.oma.be/_productsservices/coord_trans/index.php

Coordinate transformation software for France : Circé IGN

Coordinate transformation software : PROJ

Example : converting from ITRFxx (epoch XXXX) to RGF93

We have made a GNSS survey in May 2020 that we want to convert to France official referential in a projection. The GNSS referential will then be ITRF14 (epoch 2020.1), and the French referential will be RGF93 with the associated projection Lambert93. In order to do so, we need to make some referential transformation, and also some conversion between cartesian coordinates, geographic coordinates and projected coordinates.

To transform ITRF14 (epoch 2020.1) to RGF93, use this site and choose ETRF2000 (epoch 2009) as equivalent to RGF93 (see this post), for the velocities you must choose the one of the nearest IGS station (positions and velocities given at ITRF official site or Euref site). Then to transform cartesian to geographic coordinates, use Circé IGN software. And finally, use also Circé to convert to projected coordinate system.

More information

GNSS processing with open-source softwares

The main methods of GNSS processing are (1) the classical differential correction with a permanent base station, and (2) the Precise Point Positionning (PPP) method that does not require a base station and uses precises clock and ephemeris data.

The tools used can be softwares or web-services. Among software we can distinguish between commercial softwares (e.g. Trimble Business Center TBC), open-source softwares (e.g. RTKLib), and scientific softwares (Bernese, Gamit/Globk). We will introduce RTKLib and web-services.

RTKLib

RTKLib is the main open-source library and software to process GNSS data. It can do conversion, post-processing, navigation, plotting, and so on.

Note : Several versions of the software exist and they behave somewhat differently. I’ll try to describe them one by one but you’d better check carefully your results compared to another “official” reference in order to validate your process.

Versions of RTKLib :

  • Official software, with actual version 2.4.2. Does not work for me in post-processing as it is not possible to load base station Rinex data.
  • Emlid Fork here. Should be best for conversion of UBX (Emlid raw format) to Rinex. For what I can say, it worked for conversions, for post-processing of the base station, but it was not stable for post-processing of rover.
  • RTKLib_Explorer fork here. Located on the really rich blog about GNSS and RTKLib named RTKLibExplorer, this tool works well and is well documented. The setting of parameters for static processing and obtain a single position is hard (set ON for “Output Single for Sol outage” and set a value for Max Sol Std)

Complementary tools for dealing with POS files from RTKLib and CSV files from ReachView :

PPP web-services

Static GNSS correction with web-service

This processing is done with static GNSS correction using a network of permanent stations.

RGP IGN http://rgp.ign.fr/SERVICES/calcul_online.php , good for France, gives results in many coordinate systems including RGF93-Lambert93

GNSS complentary data and tools

GPS and Glonass ephemeris data, along with atmospheric data at https://kb.igs.org/hc/en-us/articles/115003935351

and also here (search more simple) https://webigs.ign.fr/gdc/fr/product/format#ephem

and also here (search simple too) https://cddis.nasa.gov/Data_and_Derived_Products/GNSS/orbit_products.html

IGS Antenna calibration file (does not contain Emlid Reach RS2) : igs14.atx

List of IGS stations : here

GPS Calendar to get the GPS week : here or at IGN site

Julian calendar to get GPS day : here or at IGN site

Offset between GPST (GPS Time) and UTC : here

More information

Filtering of DEMs

Several algorithms can be used to filter DEMs. Most simple and common are:

A less common but very recognized and efficient:

Examples of such filters on an noisy DEM :

original_dem
original_dem

Click on caption to enlarge. r : radius

Gestionnaire pour appareils Windows Mobile

La connexion directe sous Windows 10 ne fonctionne pas, il faut installer un gestionnaire PC de Windows Mobile

puis un patch correctif Windows 10

L’ensemble est disponible sur le serveur geomorpho/logiciels/terrasync/correctif_windows_mobile , ou au support technique de D3E

Pour les Juno Trimble : si l’installation des correctifs n’est pas suffisante, alors aller dans le menu du Juno Windows/Paramètres/Connexions/USBversPC et cocher (ou décocher) la case “Activer la fonctionnalité de reseau avancée”, puis attendre quelques secondes la connexion

Workflow for generating DEM from Lidar tiles

The original dataset comprises Lidar tiles with classified ground and with XYZ ASCII format of this type :

13.83535621 42.12214528 401.940 32 2
13.83534713 42.12213870 401.870 28 2
13.83533818 42.12213222 401.670 33 2
13.83532879 42.12211645 401.950 30 1
13.83533763 42.12212258 401.630 28 2
13.83534682 42.12212933 401.830 29 2
13.83535625 42.12213609 402.250 25 1

2 processing methods are proposes. The first is more straightforward, it uses bash only. The second uses ArcGIS.

Method 1 : with bash, GDAL, LAStools

All the input tiles are compressed in zip format.

1- zip extraction :

for i in *.zip; do unzip $i -d .; done

The output is  xyz files, with classification as last field. 1 for unclassified, 2 for ground (see example above)

delete zip files with :
rm *.zip

2-extract only ground points from xyz files

awk '{if ($5==2) print $0}' input > output

for batch processing, see script “script_xyz2xyzground.sh”

3-convert xyz format to las format

las2las -i input -o output

4-interpolation in grid with points2grid

for i in *.las; do BASE=`basename $i .las`;points2grid -i $i -o OUTPUT/$BASE -r 0.00000707 --input_format las --output_format arc --mean --resolution 0.00001 --fill; done

5- convert asc format to tif

for i in *.asc; do BASE=`basename $i .asc`;gdal_translate -a_srs EPSG:4326 $i ${BASE}.tif;done

6- merge all tif tiles

gdal_merge.py -o MOSAIC/Marche_lidar_full_mosaic_4326.tif -co COMPRESS=DEFLATE -co BIGTIFF=YES -co TILED=YES -a_nodata 0 *.tif

Optional step : fill holes

for i in *.tif; do BASE=`basename $i ..tif`;gdal_fillnodata.py $i ${BASE}_filled.tif -co compress=deflate -co bigtiff=yes -co tiled=yes;done

Method 2 : with ArcGIS

Toolbox of several small scripts for processing Lidar point cloud and generating DTMs

The original datasets are point clouds with ground already classified

Steps :

1 – extract only ground from the classified point cloud
This is done with awk. See description at this site https://sigeo.cerege.fr/?p=199

ex. of command for processing a file :
awk ‘{ if ( $5!=1 ) print $0 }’ input.xyz > output.xyz
to batch it in bash :
for i in *.xyz; do awk ‘{ if ( $5!=1 ) print $0 }’ input.xyz > output.xyz; done

2- Convert ASCII XYZ file to LAS
This is done with ArcGis using the txt2las tool of LAStools, but it can also be done with libLAS or PDAL
The script in the ArcGIS toolbox is used in a batch mode

3- Convert LAS files to LAS datasets
we created a small python script based on Arcpy that we use add to a Toolbox and use it in batch mode
see “convert_las2lasd.py”

4- Convert Las datasets to TIN
we created a small python script based on Arcpy that we use add to a Toolbox and use it in batch mode
see “convert_lasd2tin.py”

5- Convert Tin to raster
we created a small python script based on Arcpy that we use add to a Toolbox and use it in batch mode
see “convert_tin2raster.py”

The scripts of steps 2 to 5 are used from a ArcGIS Toolbox

Gros rasters et Compression

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

Tuilage

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 BigTIFF et les tuiles internes

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

Compression sans perte de qualité

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

Compression avec perte de qualité

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é.

Compression des pyramides raster

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 --config COMPRESS_OVERVIEW LZW dem_input.tif 2 4 8 16 32 64

Références

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

 

Traitement MicMac d’images satellites

Workflow complet pour traiter des images satellites (multi-) stéréoscopiques

avec MicMac pour le traitement de photogrammétrie, et GDAL et OTB pour les post-traitements

Génération d’images Fusionnées (Pan-Sharpen) ortho-rectifiées avec MicMac

Le process est le suivant, utilisable avec le script https://github.com/JulesFleury/JulesFleury/blob/master/Script_Decomposition_Band.sh

#0- réaliser la fusion du bundle d’origine avec otb_bundletoperfectsensor

#1-
split Pan-sharpen bands
avec otbcli_splitimage.....

#2-
copie et renommage des orientations initiales et ajustées de micmac dans le dossier Ori-XXX-adj (cf. script)

#3-
#relancer MicMac Malt pour générer les ortho sur les différents canaux, plusieurs options sont possibles

  • a- si Malt n’a pas encore tourné

#mm3d Malt Ortho "Pattern de toutes les images" "dossier_orientations" DirMEC="dossier de la MEC des panchromatiques" InMNT="Pattern des panchromatiques" DoMEC=1 DoOrtho=1 ImOrtho="Pattern des différentes bandes multispectrales" DirOF="dossier des orthophotos" EZA=1

  • b- si Malt a déjà tourné sur les panchro

#mm3d Malt Ortho "Pattern de toutes les images" "dossier_orientations" DoMEC=0 DoOrtho=1 ImOrtho="Pattern des différentes bandes multispectrales" DirOF="dossier des orthophotos pour le multispectral"
#exemple : mm3d Malt Ortho "IMG.*.TIF" Ori-RPC-d0-adj DirMEC="MEC-Malt-Zf4" DoMEC=1 ImMNT="IMG.*(00[123]).TIF" ImOrtho="IMG.*_b.*.TIF" DirOF="Ortho_Fusion" DoOrtho=1 ZoomF=4 EZA=1

#4-
#Rassembler les différents canaux des images pan-sharpen
#otbcli_ConcatenateImages ....

Lidar : Filtrage du sol et Création de MNT

Workflow pour filtrer les points du sol dans des données Lidar et créer un MNT

Ce workflow est un exemple de cas d’étude que nous avons eu à traiter, il est à adapter en fonction des données utilisées

Nous allons utiliser LibLAS, PDAL et Points2Grid

Données en entrée

Nous avons des nuages de points Lidar au format LAS, avec des coordonnées en WGS84 pour le X et Y, et en mètres pour le Z

Affichage des métadonnées (avec libLAS):

lasinfo input.las

Affichage des métadonnées (avecPDAL):

$ pdal info input.las

il faut vérifier si le système de coordonnées est définit, s’il y a un facteur d’échelle selon X, Y ou Z, et toutes autres informations

Supposons que le système de coordonnées n’est pas définit ou mal définit dans les métadonnées du LAS, cela va causer des problèmes dans la suite du traitement

Définition du système de coordonnées

Créer un fichier JSON pour le traitement PDAL, nommé par exemple “assign_scr.json” :

{
    "pipeline": [
        {
            "type" : "readers.las",
            "filename" : "input.las"
        },
        {
            "type" : "writers.las",
            "a_srs": "EPSG:4326",
            "filename" : "output.las"
        }
    ]
}

Puis lancer la commande:

pdal pipeline assign_scr.json

Maintenant la définition est correcte, mais il est préférable pour les interpolations et les traitements sur les MNT d’être en projection avec les mêmes unités en horizontal qu’en vertical, nous allons donc reprojeter en UTM

Reprojection

créer un pipeline JSON comme suit, nommé par exemple “reproject.json” :

{
  "pipeline":[
    { 
        "type":"writers.las", 
        "filename":"input.las" 
    }
    {
        "type":"filters.reprojection",
        "out_srs":"EPSG:32633"
    },
    {
      "type":"writers.las",
      "filename":"output.las"
    }
  ]
}

Puis lancer la commande:

pdal pipeline reproject.json

Il faut ensuite faire une classification sur les points

Classification du sol

créer un pipeline JSON comme suit, nommé par exemple “ground_classif.json” :

{
  "pipeline":[
    {
         "type":"filters.assign",
         "assignment":"Classification[:]=0"
    },
    {
         "type":"filters.elm"
    },
   {
         "type":"filters.outlier",
         "mean_k":8,
         "multiplier":3 
   },
   {
         "type":"filters.smrf",
         "ignore":"Classification[7:7]",
         "slope":0.2,
         "window":16,
         "threshold":0.45,
         "scalar":1.2
    },
    {
         "type":"filters.range",
         "limits":"Classification[2:2]"
     }
  ]
}

Ce pipeline est bien décrit dans les aides de PDAL, les étapes en sont les suivantes :

  1. réassigner tous les points à la classe 0 car on suppose que la classification d’origine n’est pas bonne
  2. utilisation de l’algo ELM pour assigner les points considérés comme du bruit à la classe 7
  3. faire un filtrage statistique pour les points abbérants
  4. utiliser le filtre SMRF sur les points restants pour classer les points
  5. ne récupérer que les points classés comme ground

Puis lancer la commande:

pdal pipeline ground_classif.json

on récupère donc un fichier las ne contenant que les points classés comme sol (classe 2)

il est maintenant possible d’interpoler les points pour créer un MNT

Génération du MNT

l’interpolation peut se faire avec PDAL ou avec Points2Grid (outil utilisé par OpenTopography pour générer des MNT). L’implémentation de PDAL est basée sur celle de Points2Grid, donc au final c’est équivalent je pense.

La version de points2grid que nous utilisons est celle de développement ici

points2grid -i ground.las -o dtm.tif -r 1 --output_format grid --input_format las --interpolation_mode incore --idw --resolution 1 --fill --fill_window_size 5

pour des explications sur les paramètres taper : points2grid -h

Références: