Compound CS transformation

[Midkhat]

Hi

I am involved in a drone data processing. The photos may come in GDA94 or GDA2020 and the output has to be in MGA94 with AHD geoid height.
I have figured out how to do that in Agisoft UI with the help of the manual and it works fine, but I am having trouble automating that process in API.
This is the WKT I am using when transforming GDA2020 to MGA94:

Code: [Select]

COMPD_CS["GDA94 / MGA zone 55 + AHD",
    PROJCS["GDA94 / MGA zone 55 + AHD",
        GEOGCS["GDA94",
            DATUM["Geocentric Datum of Australia 1994",
                SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],
                TOWGS84[0.06155,-0.01087,-0.04019,0.0394924,0.0327221,0.0328979,-0.009993999999999999],AUTHORITY["EPSG","9688"]],
                PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
                UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9102"]],
            AUTHORITY["EPSG","4283"]
            ],
            PROJECTION["Transverse_Mercator",AUTHORITY["EPSG","9807"]],
            PARAMETER["latitude_of_origin",0],
            PARAMETER["central_meridian",147],
            PARAMETER["scale_factor",0.9996],
            PARAMETER["false_easting",500000],
            PARAMETER["false_northing",10000000],
            UNIT["metre",1,AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","28355"]
    ],
    VERT_CS["AHD height",
        VERT_DATUM["Australian Height Datum",2005,AUTHORITY["EPSG","5111"]],
        UNIT["metre",1,AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","5711"]
    ]
]

When I try to apply this transformation to all cameras I get an error:
RuntimeError: Vertical datum missing

I have added the geoid tiff to the "C:/Program Files/Agisoft/Metashape Pro/geoids" folder, trying to add it to the Metashape.CoordinateSystem as well. Have also tried to add it to the script folder.

This is related piece of code:

Code: [Select]

   
wgs = Metashape.CoordinateSystem("EPSG::4326")

geoid_file = "C:/Program Files/Agisoft/Metashape Pro/geoids/au_ga_AUSGeoid09_V1.01.tif"
if not os.path.exists(geoid_file):           
    print(f"file {geoid_file} doesn't exist!!!!!")
mga_file = "./mine-sites/MGA55_AHD.prj"   
with open(mga_file) as f:
    wkt = f.read()
mga = Metashape.CoordinateSystem(wkt)
mga.addGeoid(geoid_file)   

for camera in cameras:
    camera.reference.location = Metashape.CoordinateSystem.transform(camera.reference.location, wgs, mga)
chunk.crs = mga
doc.save()

Best regards

[Paulo]

Hello Midkhat,

i tried defining the mga CS according to your prj file and I transformed a marker (point 1) placed in wgs84 CS to mga without any problem. See attachment...

please check that your geoid file is correctly created in geoids folder (Tiff description must correspond to VERT_DATUM entry from prj file and reference datum must be same as in prj EPSG 4283)....see 2nd attachement

Updated with latest tif geoid format where tiff description is given by from reference datum (GDA) to height datum (AHD 5711)

Last screen shot shows transformation from GDA94 / MGA zone 55 + AHD to GDA94 and resulting geoid height compared in Global mapper...

[Midkhat]

@Paulo
Thank you for your reply, I have sorted that out then and got busy with some other stuff and missed your reply, but I have recently encountered similar issue and have found your reply to my own post through search
The issue I have now - I get the Missing Vertical Datum error when import shape as Outerboundary (from .kml) and try to export raster or point cloud:
1. Export Raster (DEM or Ortho). I get error when Z component is present in .kml. As the kml is created in third party software - it is possible that it is indeed not right value. I don't really need Z in this shape, I just remove Z and this sorts it - DEM and Ortho come out clipped.
2. Export Points. In this case removing Z doesn't help as I presume it is required for point cloud clipping. At the same time I can open the project later and export clipped cloud from GUI...

Would you be able to advise or may be there is a work around for point cloud clipping?

PS Thanks again, I really appreciate your help and have learnt some new stuff from the prnscreens you shared (and other your comments in this forum)

[Paulo]

Hi MidKhat,

it is difficult to help without more info. What is the defined Compound CS you want to import your kml to? Do you have an example of typical kml you are using?