Python - set custom coordinate system

[MichaelK123]

Hello,

I am using Python for an automatic workflow. I am defining a coordinate system with

Code: [Select]

new_crs = Metashape.CoordinateSystem("EPSG::25832")
for camera in chunk.cameras:
    camera.reference.location = new_crs.project(chunk.crs.unproject(camera.reference.location))
chunk.crs = new_crs
chunk.updateTransform()
doc.save ()
which works perfectly fine.
I defined a custom coordinate system in Metashape GUI by adding a vertical cs (tiff) to the EPSG25832 (custom name: "ETRS89 / UTM zone 32N+DHHN2016". How can I set this coordinate system in Python?

Thank you in advance!
Michael

[MichaelK123]

So I did a lot of research in the Metashape forum. As suggested I tried to generate a prj.-file and a new tif-file using this manual :
https://agisoft.freshdesk.com/support/solutions/articles/31000148332-how-to-use-height-above-geoid-for-the-coordinate-system
The problem is that the newly generated tif doesn't show the "VERT-DATUM[".......]" entry in the properties (see attachments), which should be there according to this post:
https://www.agisoft.com/forum/index.php?topic=13107.0
So it seems to me that the tif wasn`t created accordingly, although the transformation with the newly generated coordinate system works fine in the Metashape GUI. If I try to change the properties of the tif-file manually, Metashape (in python) tells me that it can't load the tif-file. If I leave it without manually changing the VERT-DATUM-entry, it tells me that the vertical datum is missing.

Since this year the undulation grid for Germany (GCG2016) is available open source:
https://daten.gdz.bkg.bund.de/produkte/sonstige/quasigeoid/aktuell/quasigeoid.geo89.de.zip
Maybe someone is able to help me out out with properly creating the tif/prj-generation, or even to upload this file to the official Metashape download (https://www.agisoft.com/downloads/geoids/)?

Best regards
Michael

[Paulo]

Halo MichaelK123,

it seems the Tif provided from your link has a nodata value of  -1.#QNAN00000, which is weird. Normally it is something like 999999 as from extract from ascii GCG2016.txt file:

Code: [Select]

55.979166667  3.25625 43.2422
55.979166667  3.26875 43.2332
55.979166667  3.28125 43.2241
55.979166667  3.29375 43.2146
55.979166667  3.30625 43.2047
55.979166667  3.31875 43.1948
55.979166667  3.33125 43.1850
55.979166667  3.34375 43.1750
55.979166667  3.35625 43.1654
55.979166667  3.36875 43.1560
55.979166667  3.38125 43.1465
55.979166667  3.39375 43.1376
55.979166667  3.40625 43.1281
55.979166667  3.41875 43.1188
55.979166667  3.43125 43.1088
55.979166667  3.44375 43.0990
55.979166667  3.45625 999999 # nodata
So I would recommend you create your Metashape geoid file using the ascii format as input instead of tiff....

It should work,

[MichaelK123]

Hi Paul,

Thank you for your quick response.
As suggested I tried creating the file with the ascii format as input, still the same problem. I attached screenshots (00:settings-ASC, 01: results-ASC).
I have noticed that I have to manually select a Vertical CS when I use the ASC-format as input, but when I use the ggf-Format, Metashape automatically selects "geoid" as Vertical CS.
But even than I have the same/comparable results (screenshots 02:settings-ggf, 03: results-ggf).
Best regards
Michael

[Paulo]

Hi Michael,

I really do not know why you have still problem importing ascii file (lat/lon/height format). I redid with ascii txt file from your link and created the Compound CS with DHHN2016 height with no problem as seen from screen shot. Where I placed a point near Leipzig and calculated its geoid height and compare with same point placed on GCG2016.tif file in Global mapper. The geoid value is same...

[Paulo]

Here is the link to GCG2016.tif geoid file....https://drive.google.com/file/d/1hl4cPXq1f0rJEKir3dFV5ZXMEIzvvk8I/view?usp=drive_link

 Please copy to your Agisoft geoids folder....

[MichaelK123]

Hi Paul,

First of all, thank you very much for helping me out and creating/uploading the GCG2016-file.
I downloaded it and tried to implement it into the process. Unfortunately it is still not working...

Here is what I did:

I created a .prj-File with the downloaded GCG2016.tif file (I called it GCG2016_Paul.tif). The prj-code looks like this and I defined it as "coordinateSystem "-variable (for Metashape-Python-script):

Code: [Select]

COMPD_CS["ETRS89 / UTM zone 32N+GCG2016_Paul",PROJCS["ETRS89 / UTM zone 32N",GEOGCS["ETRS89",DATUM["European Terrestrial Reference System 1989 ensemble",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9102"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","25832"]],VERT_CS["DHHN2016",VERT_DATUM["ETRS89 (EPSG::4937) to Deutsches Haupthoehennetz 2016",2005],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]]
In Metashape-Python-script I defined the geoid path:

Code: [Select]

geoidPath = cwd + "/supportFiles/geoids/GCG2016_Paul.tif"
print("geoidPath: "+geoidPath)
and the coordinate system:

Code: [Select]

new_crs = Metashape.CoordinateSystem(coordinateSystem)
new_crs.addGeoid(geoidPath)
chunk.crs = new_crs
chunk.updateTransform()
for camera in chunk.cameras:
    camera.reference.location = new_crs.project(chunk.crs.unproject(camera.reference.location))
doc.save ()
By running the script I get the error message

Code: [Select]

Traceback (most recent call last):
  File "workflow_align2_depthMaps2_fullWorkflow.py", line 132, in <module>
    chunk.updateTransform()
RuntimeError: Vertical datum missing
I changed the name of the tiff-file (to make sure Metashape is working with the right file) and I get the error message

Code: [Select]

  File "workflow_align2_depthMaps2_fullWorkflow.py", line 130, in <module>
    new_crs.addGeoid(geoidPath)
OSError: Can't open file: No such file or directory (2): /mnt/Daten/Projekte/06_Aischtal-Staedtler/00_Datenaufnahme_2022-09-01/M17_K_h110_v4_60grad_OH60_OV60_nSP/CD/run_Flug7_full-ini/supportFiles/geoids/GCG2016_Paul.tif
I really don't get what I am doing wrong...

I am running the script on a (external) linux cluster without GUI, which is why I am not able to create the file on the cluster itself, but I thought by defining the geoid path, Metashape should know where to get the files from, right?

Best regards
Michael

[Alexey Pasumansky]

Hello Michael,

Does it make any difference, if you initialize the geoid by adding the following line before defining the coordinate system that uses the related geoid?

Code: [Select]

Metashape.CoordinateSystem.addGeoid(geoidPath)

[MichaelK123]

Hello Alexey,

unfortunately it doesn't change anything...

My code:

Code: [Select]

print('setting Geoid path:' + geoidPath)
Metashape.CoordinateSystem.addGeoid(geoidPath)
new_crs = Metashape.CoordinateSystem(coordinateSystem) #Nur ein Beispiel
#new_crs.addGeoid(geoidPath)
chunk.crs = new_crs
chunk.updateTransform()
for camera in chunk.cameras:
    camera.reference.location = new_crs.project(chunk.crs.unproject(camera.reference.location))
doc.save ()
Result:

Code: [Select]

setting Geoid path:/mnt/Daten/Projekte/06_Aischtal-Staedtler/00_Datenaufnahme_2022-09-01/M17_K_h110_v4_60grad_OH60_OV60_nSP/CD/run_Flug7_full-ini/supportFiles/geoids/GCG2016_Paul.tif
Traceback (most recent call last):
  File "workflow_align2_depthMaps2_fullWorkflow.py", line 134, in <module>
    chunk.updateTransform()
RuntimeError: Vertical datum missing
I noticed that the downloaded tif-file doesn't show the "VERT_DATUM"-entry, which is supposed to be shown according to this post: https://www.agisoft.com/forum/index.php?topic=13107.0
Could this be a problem?

Best regards
Michael

[Alexey Pasumansky]

I have transformed Paul's file via Add Geoid command in GUI, while working with compound coordinate system. Can you please check, if it works, if you replace the old file in geoids folder with the following one:
https://drive.google.com/file/d/1YHCdBDmk0HiX0U1w0YXqAEBQ7zx-wdsK/view?usp=drive_link

[MichaelK123]

Hi Alexey,

Thank you for sending the file.
I created a new .prj-file using your tiff. I moved your tiff file to the geoid folder on the linux cluster and adapted the geoid path. At the beginning it said "vertical datum out of range". The same happened when I tried to convert the coordinate system in the gui.
I reoredered some lines in the script with the result that Metashape created a new project file, importing the images and setting the coordinate system to the custom one (with the right heights!). So this is kind of a success!
Here the adapted script:

Code: [Select]

#Koordinatensystem
print('setting Geoid path:' + geoidPath)
chunk = doc.addChunk()
Metashape.CoordinateSystem.addGeoid(geoidPath)
new_crs = Metashape.CoordinateSystem(coordinateSystem)
#new_crs.addGeoid(geoidPath)
chunk.crs = new_crs

# STEP 1 - Load Images
chunk.addPhotos(photos)
doc.save()

chunk.updateTransform()
for camera in chunk.cameras:
    camera.reference.location = new_crs.project(chunk.crs.unproject(camera.reference.location))
doc.save ()
However, the calculation stops, showing the following error:

Code: [Select]

Traceback (most recent call last):
  File "workflow_align2_depthMaps2_fullWorkflow.py", line 183, in <module>
    chunk.alignCameras(adaptive_fitting = True)
RuntimeError: Infinite values in transform are not allowed
line 183 in the python script:

Code: [Select]

chunk.alignCameras(adaptive_fitting = True)
I have the feeling that there is still some problem with converting the coordinate system of the imported images (WGS84) to the custom one.
I have also tried importing the images into Metashape in GUI and convert the coordinate system to the costum one and it worked as well. So basically it can't be that wrong...

Best regards
Michael

[Alexey Pasumansky]

Hello Michael,

"Infinite values in transform are not allowed" error may indicate that there is a conflict between coordinate system selected for chunk/cameras/markers and reference coordinate values input for cameras/markers.
You can try to save the project prior to alignCameras step and then open it in GUI to see, if there is the problem with the coordinates.

[MichaelK123]

Okay, I saved the file and double checked in GUI. It seems like the images were imported with coordinates UTM32 and DHHN2016 (which is what I was aiming for), but the project coordinate system is still set to WGS84. I manually changed it to my custom coordinate system and tried to align photos. This worked!
But now I am confused. I thought I already set the right coordinate system with

Code: [Select]

Metashape.CoordinateSystem.addGeoid(geoidPath)
new_crs = Metashape.CoordinateSystem(coordinateSystem) (with the "coordinateSystem" variable referring to the input of the prj.-file)
do I have to set it again once I added the photos?

[Alexey Pasumansky]

Hello MichaelK123,

I think you need to put chunk.crs assignment lint after chunk.addPhotos() operation. If the source photos have build-in GPS tags in EXIF, the coordinates will be loaded automatically and the coordinate system will be switched to WGS84.

[MichaelK123]

You are right!
I reordered the commands/lines and now it works.

First I add the Chunk, then the photos, transform the coordinates to my custom crs, and at the end I define the chunk crs as my custom crs. Here the described part of my code:

Code: [Select]

print('setting Geoid path:' + geoidPath)
chunk = doc.addChunk()
chunk.addPhotos(photos)
Metashape.CoordinateSystem.addGeoid(geoidPath)
new_crs = Metashape.CoordinateSystem(coordinateSystem)
#new_crs.addGeoid(geoidPath)

chunk.updateTransform()
for camera in chunk.cameras:
    camera.reference.location = new_crs.project(chunk.crs.unproject(camera.reference.location))
doc.save ()

chunk.crs = new_crs

doc.save()
Thank you Alexey and Paulo for helping me out!

Best regards
Michael