Define 3 points coordinates from dip/azimuth

In case you have the dip angle and the azimuth angle (i.e. strike angle) of a plane, here is how to calculate the XYZ coordinates of 3 points of this plane.

Strike (i.e. azimuth) and dip

Note that below we assume that dip value is measured from horizontal plane. If dip is towards north, it must be set positive, and if it is towards south, set to negative.

Also, azimuth is measured positively in clockwise direction starting from North, i.e East is 90°, South is 180°…

In order to calculate coordinates of points belonging to this plane, let’s start from a point O (X0, Y0, Z0) as the origin.

Considering a distance d, a dip angle dip, and an azimuth (i.e. strike) angle az, we can use the spherical coordinates calculations to obtain the coordinates of a point A :

dXA = d*cos(az)*cos(dip)

dYA=d*sin(az)*cos(dip)

dZA=d*sin(dip)

XA=X0+dXA

YA=Y0+dYA

ZA=Z0+dZA

All angles must be set into radians.

And for a second point B on the horizontal direction (strike direction), we can do :

dXB=d*sin(az)

dYB=d*cos(az)

dZB=0

XB=X0+dXB

YB=Y0+dYB

ZB=Z0

Download a Calc file to do this here :

Usefull links to read also:

https://beckassets.blob.core.windows.net/product/readingsample/176217/9783540310549_excerpt_001.pdf

https://en.wikipedia.org/wiki/Spherical_coordinate_system#Cartesian_coordinates

https://gis.stackexchange.com/questions/13484/how-to-convert-distance-azimuth-dip-to-xyz

Batch download tiles

Imagine you need to download several tiles. Instead of interactive clicking you could save the list of tile names in a file and automatically download them with wget.

For example, to download an area of LIDAR HD IGN, here is how to do:

1- download the list of all tiles (“tableau d’assemblage”) in shapefile format

2- open it on QGis, then draw a new shapefile containing your area. Select all tiles intersecting this area and export the attribute table of selection in CSV format. You obtain a file like this:

nom_pkk,url_telech
LHD_FXX_0945_6419_PTS_O_LAMB93_IGN69.copc.laz,https://storage.sbg.cloud.ovh.net/v1/AUTH_63234f509d6048bca3c9fd7928720ca1/ppk-lidar/QN/LHD_FXX_0945_6419_PTS_O_LAMB93_IGN69.copc.laz
LHD_FXX_0953_6419_PTS_O_LAMB93_IGN69.copc.laz,https://storage.sbg.cloud.ovh.net/v1/AUTH_63234f509d6048bca3c9fd7928720ca1/ppk-lidar/QN/LHD_FXX_0953_6419_PTS_O_LAMB93_IGN69.copc.laz
LHD_FXX_0959_6422_PTS_O_LAMB93_IGN69.copc.laz,https://storage.sbg.cloud.ovh.net/v1/AUTH_63234f509d6048bca3c9fd7928720ca1/ppk-lidar/QN/LHD_FXX_0959_6422_PTS_O_LAMB93_IGN69.copc.laz

3- extract only the URLs from this file with awk

awk -F, '$2~/https/ {OFS=",";print $2}' input.csv > urls.csv

4- download all urls with wget

wget -i urls.csv

This is it! Maybe one should add a verification step to be sure that all files are properly downloaded.

EarthAnalytics Python installation

voir https://www.earthdatascience.org/workshops/setup-earth-analytics-python/

il faut installer git, conda (préferer miniconda) et bash

ensuite on peut créer l’environnement earthanalytics en utilisant le fichier yaml et la méthode décrite ici .

pour lister les environnements disponibles taper

conda env list

pour activer l’environnement taper

conda activate <nom_environnement>

Plus de détails sur la CheatSheet de Conda ici

Ensuite il faut choisir un éditeur de texte. Dans mon cas, vu que j’ai installé miniconda et l’environnement virtuel sous Linux en WSL, j’utilise Spyder depuis Windows en définissant le kernel, cf. install Spyder on WSL2 Linux , et notamment cette commande depuis le Shell

python -m spyder_kernels.console --matplotlib="inline" --ip=127.0.0.1 -f=~/remotemachine.json &

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/

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

examples

  • 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

  • from Lambert93 planimetric and NGF-IGN69 (with RAF20) to WGS84 ellipsoidal heights

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

  • 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