Ellipsoidal to orthometric altitudes for geoid model

The geoid models have sometimes a hard to decode format, such as IGN .mnt files https://geodesie.ign.fr/index.php?page=grilles

For now, I’ve just found third parties software that can do the format conversion, even if this could be doable with a bash or python script.

Such files can be converted to ASCII XYZ format with the demo version of Hydromagic, an hydrographic surveying software. Here is how to do :

https://www.eye4software.com/hydromagic/documentation/manual/utilities/geoid-file-conversion/

CRX to RNX

CRX is a Hatanaka compression for Rinex file used by many permanent networks. Y. Hatanaka of GSI wrote softwares for compress/uncompress, see https://terras.gsi.go.jp/ja/crx2rnx.html.

To use the software (example of CRX compression to uncompressed RNX format):

  • open a shell (e.g. Powershell)
  • go to the folder containing the executable
  • type “CRX2RNX.exe input.crx”
  • you get as output a RNX file that you can rename to .yyo
  • for help type “CRX2RNX.exe -h”

Use Reach M2 with Phantom 4 Pro

dataset :

photos with timestamps

UBX from M2 and from base station

Processing:

1- Post-Process base station and get corrected coordinates

2- with Emlid Studio (or RTKLib) calculate events (triggers) positions from UBX

as output one gets an …._events.pos file that gives corrected coordinates for each event together with a timestamp

3- extract timestamps from exif of pictures

use exiftool with this command

exiftool -n -csv -filemodifydate -api TimeZone=UTC *.JPG > metadata.csv

the output file “metadata.csv” gives the name of the picture and the timestamp in UTC format

4- Merge …_events.pos and pictures events

For now I’ve got this ugly Python script:

#! /usr/bin/env python

"""
Update Emlid Reach Survey points with PPK position output from RTKLIB
David Shean
dshean@gmail.com
Edited to fix Pandas datetime/Timestamp tz issues, and a few key changes likely based on Emlid updates
"""

import os
import argparse
import numpy as np
import pandas as pd

def getparser():
    parser = argparse.ArgumentParser(description='Update Emlid Reach Survey points with \
            PPK positions from RTKLIB')
    parser.add_argument('survey_pts_csv_fn', type=str, help='Survey point csv filename')
    parser.add_argument('ppk_pos_fn', type=str, help='PPK pos filename')
    return parser


def main():
    parser = getparser()
    args = parser.parse_args()

    survey_pts_csv_fn = args.survey_pts_csv_fn
    ppk_pos_fn = args.ppk_pos_fn

    print('Loading: %s' % survey_pts_csv_fn)
    survey_pts = pd.read_csv(survey_pts_csv_fn, parse_dates=[1], header=0)
    survey_pts['date']=pd.to_datetime(survey_pts['FileModifyDate'],format="%Y:%m:%d %H:%M:%S+00:00")
    survey_pts.sort_values('date', inplace=True)
    survey_pts.index=survey_pts['date']
    print(survey_pts.dtypes)
    print(survey_pts)
    header = 'Date UTC latitude(deg) longitude(deg)  height(m)   Q  ns   sdn(m)   sde(m)   sdu(m)  sdne(m)  sdeu(m)  sdun(m) age(s)  ratio'
    print('Loading: %s' % ppk_pos_fn)
    ppk_pos = pd.read_csv(ppk_pos_fn, comment='%', delim_whitespace=True, names=header.split(), parse_dates=[[0,1]])
    ppk_pos['date']=pd.to_datetime(ppk_pos['Date_UTC'])
    ppk_pos.index=ppk_pos['Date_UTC']
    print(ppk_pos.dtypes)
    print(ppk_pos)

    # Applying merge_asof on data and store it
    # in a variable
    merged_dataframe = pd.merge_asof(ppk_pos, survey_pts, right_index=True,left_index=True,direction='nearest',tolerance=pd.Timedelta("1s"))

    print(merged_dataframe)

    #Write out new file
    out_fn = os.path.splitext(survey_pts_csv_fn)[0]+'_merged.csv'
    print("Writing out: %s" % out_fn)
    merged_dataframe.to_csv(out_fn)    

if __name__ == "__main__":
    main()

5- In Metashape, import the coordinates of the cameras from the new file

Convert raster from ellispoidal height to orthometric altitude

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/

using gdalwarp to convert from UTM33N-WGS84 ellipsoidal height to UTM33N-EGM96 geoid height

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

The definition of the projection is the PROJ4 line that can be found on https://spatialreference.org/

Inversely, to convert from altitude to ellipsoidal height use such command (change srs and geoid grid if needed):

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

Trame NMEA GeoXT

Pour activer la sortie des trames NMEA sur le GeoXT Trimble :

  • installer l’adaptateur RS232
  • Brancher un cable RS232 (attention certains cables ne fonctionnent pas)
  • Menu du GeoXT / Paramètres / Connexion / GPS Connector : cela lance la sortie des trames NMEA sur le port COM.
  • Les paramètres de sortie sont accessibles

Metashape script

API Python voir https://www.agisoft.com/pdf/metashape_python_api_1_5_0.pdf

FAQ scripting : https://agisoft.freshdesk.com/support/solutions/folders/31000114192

Liste des scripts: https://github.com/agisoft-llc/metashape-scripts

exemple d’utilisation de la console Python dans Metashape, décocher les orientations de caméras :

doc=Metashape.app.document
chunk=doc.chunk
for camera in chunk.cameras:
   camera.reference.rotation_enabled=0

GPS time

GPS time system (recording the absolute number of seconds since the start of the GPS epoch at midnight on January 6th 1980)

Sometimes GPS time is given with an offset (e.g. 1000000000), in this case the 0 value corresponds to 2011/09/14 at 01:46:25 UTC.

Note also that POSIX (Unix) time is given since Jan 1st, 1970, 00:00:00 UTC.

convert to date

835734.87 6250812.50 -3.04 20 0 57560982
835734.88 6250296.85 -8.16 30 0 46351553
835734.88 6250476.48 -6.97 30 0 57559145

example of gps dataset

The gps epoch is in the sixth field, and the values are offset by 1 billion. To convert to date one can use Python.

>>> import datetime
>>> datetime.date.fromtimestamp(1315964785) #posix time of the GPS epoch including the 1 billion offset
datetime.date(2011, 9, 14)
>>> datetime.date.fromtimestamp(1315964782+57560982) # gps time of the dataset + the posix time
datetime.date(2013, 7, 11)

Tachéo Leica TC805

Téléchargement données instrument

Windows 7

On peut télécharger les données acquises par le tachéomètre Lecia TCR 805 avec le logiciel Leica Geo Office. Ce logiciel ne s’installe que sur Windows 7. L’installation sur Windows 10 n’aboutit pas. Le logiciel se trouve ici: \\geomorpho\Logiciel\leica\Leica Geo Office Tools_fr.zip
Le téléchargement des données est très simple avec un câble série(pas d’installation de driver requise). Elle peut poser des problèmes avec un cable USB.

Pour le cable USB, le driver s’installe tout seul (validé sur Windows 7). Ensuite utiliser l “Data Exchange Manager”, il faut sur le port COM correspondant configurer les paramètres avec l’instrument TPS800 et le taux de transfert 19200 baud. Les autres paramètres peuvent être laissés par défaut.

Ensuite dans la fenêtre de transfert des données, il suffit de faire glisser les fichiers, et choisir le format IDX.

Windows 10

Avec le cable USB, le driver s’installe tout seul (validé sur Windows 7 et aussi sous Windows 10). Ensuite utiliser en tant qu’Administrateur le logiciel “Data Exchange Manager”, il faut sur le port COM correspondant configurer les paramètres avec l’instrument TPS800 et le taux de transfert 19200 baud. Les autres paramètres peuvent être laissés par défaut.

Ensuite dans la fenêtre de transfert des données, il suffit de faire glisser les fichiers, et choisir le format IDX.

Manuels

docs dispos ici

Manuel en français plus succinct
https://www.google.fr/url?sa=t&source=web&rct=j&url=http://www.grr.ulaval.ca/gae_3005/Labos/Docs/CH_Releves.pdf&ved=2ahUKEwi56OrhyPTnAhUU4OAKHbenBbcQFjASegQIBhAB&usg=AOvVaw0o16NjveUM588rvHZfshAJ

Manuel complet
https://www.google.fr/url?sa=t&source=web&rct=j&url=https://tmackinnon.com/2005/PDF/Leica_tc605_tc805_tc905_uesr_manual.pdf&ved=2ahUKEwiX2fOHyPTnAhVHD2MBHYBGAMgQFjAAegQIAhAB&usg=AOvVaw3p6oIfs2I_jrw8G1qsnLIA

Crop RPC

Using MicMac tools :

example:

mm3d SateLib CropRPC Ori-RPC-d0/GB-Orientation-IMG_PHR1A_P_201902190719128_SEN_3788788101-001.TIF.xml Ori-RPC-d0/GB.* Cropped Org=[3000,4000] Sz=[10500,8500]

parameters:

  1. image to use for the definition of crop zone (Org and Sz correspond to this image)
  2. pattern of orientation files for images to be cropped
  3. name of folder to store cropped RPC
  4. Org : Origin of the box to crop
  5. Sz : Size of the box to crop