Image footprint coverage #850
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,155 @@ | ||
| import math | ||
|
|
||
| import pyproj | ||
| from shapely.affinity import rotate | ||
| from shapely.geometry import Polygon, mapping, shape | ||
| from shapely.ops import transform, unary_union | ||
|
|
||
|
|
||
| # Convert lon/lat to meters for footprint and area math. | ||
| # A degree is not a fixed ground distance, so we avoid area math in EPSG:4326. | ||
| projector = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) | ||
| inverse_projector = pyproj.Transformer.from_crs( | ||
| "EPSG:3857", "EPSG:4326", always_xy=True | ||
| ) | ||
|
|
||
| # Fallbacks when EXIF/project metadata is missing. | ||
| DEFAULT_DIAGONAL_FOV_DEG = 82.1 | ||
| DEFAULT_IMAGE_WIDTH = 4000 | ||
| DEFAULT_IMAGE_HEIGHT = 3000 | ||
|
|
||
|
|
||
| def _as_float(value) -> float | None: | ||
| # EXIF values can be messey as strings, blanks, or missing. | ||
| # This function converts values into numbers or None. | ||
| if value is None: | ||
| return None | ||
|
|
||
| try: | ||
| return float(value) | ||
| except (TypeError, ValueError): | ||
| return None | ||
|
|
||
|
|
||
| def _footprint_size_meters( | ||
| gsd_cm_px: float | None, | ||
| altitude_m: float | None, | ||
| image_width: float | None, | ||
| image_height: float | None, | ||
| ) -> tuple[float, float] | None: | ||
| """Estimate the ground width/height covered by one image.""" | ||
| # Prefer GSD: it tells us directly how much ground each pixel covers. | ||
| gsd_cm_px = _as_float(gsd_cm_px) | ||
| altitude_m = _as_float(altitude_m) | ||
| image_width = _as_float(image_width) or DEFAULT_IMAGE_WIDTH | ||
| image_height = _as_float(image_height) or DEFAULT_IMAGE_HEIGHT | ||
|
|
||
| if gsd_cm_px and gsd_cm_px > 0: | ||
| # Example: 2 cm/px * 4000 px / 100 = 80m on the ground. | ||
| return ( | ||
| gsd_cm_px * image_width / 100, | ||
| gsd_cm_px * image_height / 100, | ||
| ) | ||
|
|
||
| # If no GSD: altitude + camera FOV gives an approximate ground rectangle. | ||
| if altitude_m and altitude_m > 0: | ||
| # Keep the rectangle shaped like the image, usually 4:3. | ||
| aspect_ratio = image_width / image_height | ||
| diagonal_m = ( | ||
| 2 * altitude_m * math.tan(math.radians(DEFAULT_DIAGONAL_FOV_DEG) / 2) | ||
| ) | ||
| # Convert the estimated diagonal into width and height. | ||
| height_m = diagonal_m / math.sqrt(1 + aspect_ratio**2) | ||
| width_m = aspect_ratio * height_m | ||
| return width_m, height_m | ||
|
|
||
| return None | ||
|
|
||
|
|
||
| def image_footprint_polygon(image: dict) -> Polygon | None: | ||
| """Build one image footprint as a yaw-rotated rectangle in meters.""" | ||
| # No GPS point means no map footprint. | ||
| if not image.get("location"): | ||
| return None | ||
|
|
||
| # Calculate how big this image is on the ground. | ||
| # GSD preferred, altitude/FOV is the second choice. | ||
| footprint_size = _footprint_size_meters( | ||
| image.get("gsd_cm_px"), | ||
| image.get("altitude_m"), | ||
| image.get("image_width"), | ||
| image.get("image_height"), | ||
| ) | ||
| if not footprint_size: | ||
| return None | ||
|
|
||
| # Build in meters first, then convert back to lon/lat only for display. | ||
| width_m, height_m = footprint_size | ||
| # Convert image GPS point into meter coordinates | ||
| center = transform(projector.transform, shape(image["location"])) | ||
| x, y = center.x, center.y | ||
| half_width = width_m / 2 | ||
| half_height = height_m / 2 | ||
|
|
||
| # Draw rectangle corners (not rotated) | ||
| footprint = Polygon( | ||
| [ | ||
| (x - half_width, y - half_height), | ||
| (x + half_width, y - half_height), | ||
| (x + half_width, y + half_height), | ||
| (x - half_width, y + half_height), | ||
| (x - half_width, y - half_height), | ||
| ] | ||
| ) | ||
|
|
||
| # Rotate it using drone/camera heading when EXIF has yaw. | ||
| yaw_deg = _as_float(image.get("yaw_deg")) | ||
| if yaw_deg is not None: | ||
| # Shapely rotates counter-clockwise from east; yaw is clockwise from north. | ||
| footprint = rotate(footprint, 90 - yaw_deg, origin=center) | ||
|
|
||
| return footprint | ||
|
|
||
|
|
||
| def image_footprints_feature_collection(images: list[dict]) -> dict: | ||
| """Return map-ready GeoJSON outlines for all image footprints.""" | ||
| # Frontend map layers need GeoJSON in lon/lat. | ||
| features = [] | ||
| for image in images: | ||
| footprint = image_footprint_polygon(image) | ||
|
Member
There was a problem hiding this comment. **any way we could consolidate these two footprint calcs, so we only have to do it once per image? |
||
| if footprint is None: | ||
| continue | ||
|
|
||
| # Convert meters back to lon/lat. MapLibre expects GeoJSON coordinates in that. | ||
| footprint_wgs84 = transform(inverse_projector.transform, footprint) | ||
| features.append( | ||
| { | ||
| "type": "Feature", | ||
| "geometry": mapping(footprint_wgs84), | ||
| "properties": {"image_id": str(image["id"])}, | ||
| } | ||
| ) | ||
|
|
||
| return {"type": "FeatureCollection", "features": features} | ||
|
|
||
|
|
||
| def coverage_percentage_from_footprints( | ||
|
Member
There was a problem hiding this comment. This is great 👍 We calculate coverage in this endpoint, but also for the entire project in It would be great to refactor the usage there to also use this updated code, so the project level coverage is calculated correctly (merging this as-is would cause two different values between the two implementations. The old implementation can be removed / cleaned up) |
||
| coverage_geometry: dict, | ||
| images: list[dict], | ||
| ) -> float: | ||
| """Calculate covered area percentage from unioned rectangular footprints.""" | ||
| # Merge rectangles, clip to the task, then divide by task area. | ||
| target_m = transform(projector.transform, shape(coverage_geometry)) | ||
| footprints = [ | ||
| footprint | ||
| for footprint in (image_footprint_polygon(image) for image in images) | ||
| if footprint is not None | ||
| ] | ||
|
|
||
| if not footprints or target_m.is_empty or target_m.area <= 0: | ||
| return 0.0 | ||
|
|
||
| # unary_union counts overlapping photo footprints only once | ||
| # photos outside the task area don't count | ||
| covered_m = unary_union(footprints).intersection(target_m) | ||
| return min(100.0, (covered_m.area / target_m.area) * 100) | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This will be an EPSG:3857 (web mercator) coordination, but the footprint calculated below is true ground meters, so the calculations as we get closer to the equator will be increasingly incorrect.
I would propose a simple solution is to correct by the
1/cos(latitude)adjustment that web mercator uses: