Forum

Author Topic: Extract conents of report using Python  (Read 1402 times)

ankit108

  • Newbie
  • *
  • Posts: 17
    • View Profile
Extract conents of report using Python
« on: February 18, 2023, 08:57:32 AM »
Hello all. I am trying to fetch the contents of the results using Python, that are there when we do "Generate Report" in Agisoft.
Usually on the 2nd page of report, we have "Survey Data". I want to get the contents of that. More specifically, I want to fetch the value corresponding to "Coverage Area" and Resolution.
Basically I am making a csv file, where I am putting the values of total Cameras, Total RMSE Projection Error, Total Marker Error, Resolution and Coverage Area.
So far I have succeeded in getting total no. of cameras (easy), RMSE for cameras. Total error of marker, I am still trying, with little success, but for Resolution and "Coverage Area", I have no clue.
Kindly someone help me regarding that.. Please..
« Last Edit: February 19, 2023, 10:12:51 AM by ankit108 »

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: Extract conents of report using Python
« Reply #1 on: February 20, 2023, 03:44:09 PM »
Hello ankit108,

The ground sampling resolution can be calculated as flight altitude in meters divided by the focal length in pixels. The latter can be calculated as focal length in mm divided by the pixel size on the sensor.

To estimate the flight altitude value you can perform the following operation:
at first estimate the average distance from the camera center to the tie points that have valid 2D projections on the given photo and then takes the median value for all the aligned cameras in the project.
Then the ground resolution value (as for the report file) can be estimated individually for each camera as the flying altitude (like previously described) multiplied by the angular resolution of the pixel. Then the value is averaged for all cameras.

 
The following script calculates GSD and flight height according to the provided description (script for 2.0 version):
Code: [Select]
import Metashape, statistics, math

def percentile(data, perc: int):
size = len(data)
return sorted(data)[int(math.ceil((size * perc) / 100)) - 1]

def estimate_gsd_height(chunk):

cameras = [camera for camera in chunk.cameras if camera.transform and camera.type == Metashape.Camera.Type.Regular]
if not len(cameras):
print("No aligned cameras, script aborted")
return 0

point_cloud = chunk.tie_points
points = point_cloud.points
npoints = len(points)
projections = point_cloud.projections
T = chunk.transform.matrix

cameras = [camera for camera in list(cameras) if (camera in projections)]

point_ids = [-1] * len(point_cloud.tracks)
for point_id in range(0, npoints):
point_ids[points[point_id].track_id] = point_id

dist = list()
flight_height = list()

for camera in cameras:
sensor = camera.sensor
f = sensor.calibration.f
dist_camera = list()
altitude = list()

for proj in projections[camera]:
track_id = proj.track_id
point_id = point_ids[track_id]
if point_id < 0:
continue
point = points[point_id]
if not point.valid:
continue

coord = point.coord
coord.size = 3
d = (T.mulp(coord) - T.mulp(camera.center)).norm()
dist_camera.append(d / f)
altitude.append(d)

flight_height.append(percentile(altitude, 20))
#print(camera.label, flight_height[-1])
#flight_height.append(statistics.median(altitude))
dist.append(statistics.median(dist_camera))

X = statistics.median(flight_height)
#X = sum(flight_height) / len(flight_height)
gsd = sum(dist)/len(dist)
print("GSD = " + str(gsd))
print("Flight height = " + str(X))
return gsd, X

chunk = Metashape.app.document.chunk
estimate_gsd_height(chunk)
Best regards,
Alexey Pasumansky,
Agisoft LLC

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: Extract conents of report using Python
« Reply #2 on: February 20, 2023, 03:49:57 PM »
As for the area covered, currently there's no direct way to get this value via Python, but there are some workarounds. One of our users suggested the following option, based on the generated Seamlines with zero-simplification:

Code: [Select]
# Function for measuring orthomosaic area
def get_ortho_area(chunk):

    # Build seamlines
    chunk.buildSeamlines(epsilon = 0)

    # Get shapes from seamlines group and add area
    tot_area = 0
    for shape in list(chunk.shapes):
        if shape.group.label == "Seamlines":
            tot_area = tot_area + shape.area()
            chunk.shapes.remove(shape)

    # Remove seamlines
    chunk.shapes.remove([group for group in chunk.shapes.groups if group.label == "Seamlines"])

    return int(tot_area)
Best regards,
Alexey Pasumansky,
Agisoft LLC

ankit108

  • Newbie
  • *
  • Posts: 17
    • View Profile
Re: Extract conents of report using Python
« Reply #3 on: February 20, 2023, 07:26:53 PM »
Thank you so much Alexey..
I will apply and check it in my project and will let you know.
Thanks a lot.
I have been working to solve this for last one month.

My idea was (simpler):
1. Generate footprint of images (code on github)
2. Take the union of all the shapes and then find the area of the union.

My 2nd idea was (complicated):
1. Generate image footprint
2. add area of shapes and subtract area of intersection from it.

I even took help from chatGPT, but that too didn't work.

Similarly, I am working on one more problem.
Just like image footprint, I wanted to generate footprint of Ortho also.
But this one is still unsolved.
I will be highly obliged, if you can provide any help on that as well.

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: Extract conents of report using Python
« Reply #4 on: February 20, 2023, 07:32:53 PM »
Hello ankit108,

I believe that the procedure with the seamlines would be faster and easier to implement, as multiple footprints can overlap and you will need to consider that.

And what do you mean by the footprint of Orthomosaic? The shape that follows the boundary of the used pixels?
Best regards,
Alexey Pasumansky,
Agisoft LLC

ankit108

  • Newbie
  • *
  • Posts: 17
    • View Profile
Re: Extract conents of report using Python
« Reply #5 on: February 21, 2023, 01:32:11 PM »
Yes, Exactly.. footprint that follows the boundary of Orthomosaic (visible pixels, i.e. non-nodata pixels)..
Is it doable in Metashape?? either in generic way or with some workarounds??
if yes, kindly provide some way..
It will be of a great help for me..
Thank You !!

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: Extract conents of report using Python
« Reply #6 on: February 21, 2023, 01:48:35 PM »
Hello ankit108,

At the moment there is no built-in feature to get the orthomosaic boundary from Python API, so the workaround could be a concatenation of the generated seamlines layer.
Best regards,
Alexey Pasumansky,
Agisoft LLC

ankit108

  • Newbie
  • *
  • Posts: 17
    • View Profile
Re: Extract conents of report using Python
« Reply #7 on: February 22, 2023, 09:20:06 AM »
Hello Alexey,
Thank you for your reply.

After taking a hint from your answer, just for the sake of experiment, I processed 250 Cameras, Coverage Area = 0.86 Sq. Km.
I used same code that you posted here to calculate Coverage Area, I just commented out the lines, where you are removing the shapes, i.e. Seamlines..
Then I exported the shapes. I thought it will export polygons in multiple layers, but that did not happen, so I thought, it worked.
But when I opened that in a GIS software, it had total. 29162 features.
So, theoretically, it worked, but practically it is not feasible. because I generally have 1500-2500 images per flight, which will have around 150k+ - 200k+ features in Seamlines. Running any spatial query on that will be tough and undesirable.

As you mentioned "Concatenation of Seamline Features" in your answer, can it be done in Agisoft?
Is there any way in Metashape, where in something like a for loop, we can concatenate seamlines to one-another, which in turn will make it to a big polygon with only one feature (may be).
It will become like a footprint of ortho in true sense.


If it is doable in that way, then OMG, my weeks of search will finally come to an end.
Can you please help me with that, sir?