Forum

Show Posts

This section allows you to view all posts made by this member. Note that you can only see posts made in areas you currently have access to.


Messages - Midkhat

Pages: [1]
1
Python and Java API / Exports in local coordinates
« on: October 24, 2022, 04:35:02 AM »
I am working on a project where the export is required in local grid.
The cameras data comes in GDA2020 and I have a set of GCPS in MGA2020 and corresponding in Local Grid, something like that:
MGA:
682963.000,6290736.000,751.946
688160.278,6297960.584,1000.842
680000.000,6290000.000,1000.000
LOCAL:
14586.825,15604.061,5751.946
15163.831,24486.096,6000.842
12464.381,13408.995,6000.000

I know it can be done with setting TOWGS84 Helmert 7-parameters in Agisoft, but struggle to figure out how I get these 7 parameters. Could someone direct me to how I can build a python script to define these parameters from the GCPs?

At the moment I solve that with transforming rasters and point clouds from MGA to Local with affine transformation, which generally works OK, but there are benefits of exporting in local directly from Agisoft in our project.

2
I have tried to use the Geoid which is in Applications/Metashape/../geoids folder, same processing fails at Match photos:
Code: [Select]
2022-07-24 23:38:08,579 - DEBUG - Geoid path: /Applications/MetashapePro.app/Contents/MacOS/geoids/au_ga_AUSGeoid09_V1.01.tif
2022-07-24 23:38:08,579 - DEBUG - Geoid EPSG: '4939'
2022-07-24 23:38:08,579 - DEBUG - Geoid description: 'GDA94 (EPSG::4939) to AHD height (EPSG::5711)'
2022-07-24 23:38:08,580 - DEBUG - Chunk CRS set to
2022-07-24 23:38:08,580 - DEBUG - COMPOUND WKT:

COMPD_CS["GDA94 + AHD height",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","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9102"]],AUTHORITY["EPSG","4283"]],VERT_CS["AHD height",VERT_DATUM["Australian Height Datum",2005,AUTHORITY["EPSG","5111"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","5711"]],AUTHORITY["EPSG","9464"]]

2022-07-24 23:38:08,787 - INFO - Processing step 'MatchPhotos' started.
MatchPhotos: accuracy = Medium, preselection = generic, reference, keypoint limit = 40000, keypoint limit per mpx = 1000, tiepoint limit = 4000, apply masks = 0, filter tie points = 0, filter stationary points = 1, guided matching = 1
2022-07-24 23:38:08,997 - CRITICAL - matchPhotos failed.
2022-07-24 23:38:08,997 - CRITICAL - Vertical datum missing

3
I need to set chunk.crs to 'GDA + AHD', which does exist inbuilt set of CRSs. When using API I need to be able to use geoids from some particular folder, not Metashape one, so I have exact same geoids in Metashape geoid folder and in another one.
When I set CRS in UI (either through Reference Settings menu or with console)  and run Alignment - all works fine.
When I do the same with API - set chunk CRS to 'GDA94 + AHD' and add geoid - MatchPhotos/Alignment fails with the error: "Vertical datum missing". But if I open the file after that - CRS is set to 'GDA94 + AHD' with geoid and Alignment from UI works fine again...
I have checked that the files in Metashape and custom geoids folder are identical.
I have attached two screenshots:
1. Setting CRS in UI - alignment runs fine
2. Opening document after setting CRS from standalone script and failing alignment (if align from this point in UI - works fine)
Code:
Code: [Select]
    crs = [x for x in Metashape.CoordinateSystem().listBuiltinCRS() if x.name == 'GDA94 + AHD height'][0]
    log.debug(f"CRS WKT: {crs.wkt}\n")

    #### get geoid info
    data = gdal.Open(geoid_path)
    metadata = data.GetMetadata()
    wkt = data.GetProjection()
    proj = osr.SpatialReference(wkt=wkt)
    epsg_tiff = proj.GetAttrValue('AUTHORITY', 1)
    log.debug(f"Geoid path: {geoid_path}")
    log.debug(f"Geoid WKT: {wkt}")
    log.debug(f"Geoid EPSG: {epsg_tiff}")
    log.debug(f"Geoid CRS name: {crs.name}")
    log.debug(f"Geoid description: {metadata['TIFFTAG_IMAGEDESCRIPTION']}")
    log.debug(f"CRS WKT: \n{crs.wkt}\n")

    crs.addGeoid(geoid_path)

    chunk.crs = crs
    log.info(f"Chunk <{chunk.label}> CRS is '{chunk.crs.name}'")

Response:
Code: [Select]
2022-07-23 16:59:34,506 - DEBUG - CRS WKT:
COMPD_CS["GDA94 + AHD height",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","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9102"]],AUTHORITY["EPSG","4283"]],VERT_CS["AHD height",VERT_DATUM["Australian Height Datum",2005,AUTHORITY["EPSG","5111"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","5711"]],AUTHORITY["EPSG","9464"]]

2022-07-23 16:59:34,524 - DEBUG - Geoid path: ./metashape/geoids/au_ga_AUSGeoid09_V1.01.tif
2022-07-23 16:59:34,524 - DEBUG - Geoid WKT: GEOGCRS["GDA94",DATUM["Geocentric Datum of Australia 1994",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1]],USAGE[SCOPE["Geodesy."],AREA["Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],BBOX[-60.55,93.41,-8.47,173.34]],ID["EPSG",4939]]
2022-07-23 16:59:34,524 - DEBUG - Geoid EPSG: 4939
2022-07-23 16:59:34,524 - DEBUG - Geoid CRS name: GDA94 + AHD height
2022-07-23 16:59:34,524 - DEBUG - Geoid description: GDA94 (EPSG::4939) to AHD height (EPSG::5711)
2022-07-23 16:59:34,524 - DEBUG - CRS WKT:
COMPD_CS["GDA94 + AHD height",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","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9102"]],AUTHORITY["EPSG","4283"]],VERT_CS["AHD height",VERT_DATUM["Australian Height Datum",2005,AUTHORITY["EPSG","5111"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","5711"]],AUTHORITY["EPSG","9464"]]

2022-07-23 16:59:34,525 - INFO - Chunk <API_chunk> CRS is 'GDA94 + AHD height'

2022-07-23 16:59:34,933 - CRITICAL - matchPhotos failed.
2022-07-23 16:59:34,933 - CRITICAL - Vertical datum missing



4
Python and Java API / Re: Compound CS transformation
« on: March 29, 2022, 02:14:24 AM »
@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)

5
Bug Reports / failed: error 1286
« on: January 24, 2022, 03:31:57 AM »
I get bunch of following two errors since switched to 1.8:

2022-01-24 11:13:01 OpenGL Vendor: Apple
2022-01-24 11:13:01 OpenGL Renderer: Apple M1
2022-01-24 11:13:01 OpenGL Version: 2.1 Metal - 76.3
2022-01-24 11:13:01 Maximum Texture Size: 16384
2022-01-24 11:13:01 Quad Buffered Stereo: not enabled
2022-01-24 11:13:01 ARB_vertex_buffer_object: supported
2022-01-24 11:13:01 ARB_texture_non_power_of_two: supported
2022-01-24 11:13:01 Agisoft Metashape Professional Version: 1.8.0 build 13794 (64 bit)
2022-01-24 11:13:01 Platform: Mac OS
2022-01-24 11:13:01 RAM: 16.0 GB
2022-01-24 11:13:02 LoadProject: path =
2022-01-24 11:13:02 loaded project in 0.008163 sec
2022-01-24 11:13:02 Finished processing in 0.008331 sec (exit code 1)
2022-01-24 11:13:02 renderer_qgl.cpp line 659: blitFramebuffer(0, 0, width, height, 0, 0, width, height, 0x00000100, 0x2600) failed: error 1286
2022-01-24 11:13:02 renderer_qgl.cpp line 674: glBlitFramebuffer(0, 0, width, height, 0, 0, width, height, 0x00000100, 0x2600) failed: error 1286
....
2022-01-24 11:13:05 renderer_qgl.cpp line 659: blitFramebuffer(0, 0, width, height, 0, 0, width, height, 0x00000100, 0x2600) failed: error 1286
2022-01-24 11:13:05 renderer_qgl.cpp line 674: glBlitFramebuffer(0, 0, width, height, 0, 0, width, height, 0x00000100, 0x2600) failed: error 1286

They only appear when I create a shape, or import a kml (in both cases - from GUI or API), save and reopen the project.
The errors keep streaming in even when I only move a mouse (or Pan, Zoom, click) in Model view. They stop coming in when switch to Ortho view.

6
Python and Java API / Compound CS transformation
« on: December 15, 2021, 02:00:40 AM »
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

Pages: [1]