Skip to content

Commit

Permalink
fix: Make capture area geometry valid TDE-1306 (#1163)
Browse files Browse the repository at this point in the history
Caveat: Shapely's `BaseGeometry.buffer` produces different results on
different machines, even with exactly the same Python package versions
installed. To reproduce, change `test_should_make_valid_capture_area` to
run with a buffer distance of `0.3`, and revert the change in
`merge_polygons`. At this point the test passes on my machine, but fails
on @schmidtnz's machine! We tracked this down to the line
`union_buffered.buffer(-buffer_distance, cap_style=BufferCapStyle.flat,
join_style=BufferJoinStyle.mitre)`. With the changes mentioned above,
this returns a `MultiPolygon` on my machine (basically doing the same as
`make_valid`, and therefore making the test pass), and a `Polygon` on
@schmidtnz's machine (making the `make_valid` necessary to produce a
valid polygon).

### Motivation

Ensure that published `capture-area.geojson` files are valid.

### Modifications

Use [Shapely's
`make_valid`](https://shapely.readthedocs.io/en/2.0.4/reference/shapely.make_valid.html)
before saving the geometry.

### Verification

- [x] Unit test
- [x] Verify in QGIS:
   1. `nix-shell --packages qgis --run qgis`
   3. Press <kbd>Ctrl-l</kbd>
4. Open one of the `capture-area.geojson` files created by the unit test
(for example by pausing the test at the assert)
   5. Go to Vector → Geometry Tools → Check Validity
   6. Press _Run_
   7. Verify that it prints `INVALID_COUNT: 0` and `ERROR_COUNT: 0`
  • Loading branch information
l0b0 authored Nov 7, 2024
1 parent df5a746 commit be8a039
Show file tree
Hide file tree
Showing 2 changed files with 25 additions and 1 deletion.
3 changes: 2 additions & 1 deletion scripts/stac/imagery/capture_area.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@

from linz_logger import get_log
from shapely import BufferCapStyle, BufferJoinStyle, to_geojson, union_all, wkt
from shapely.constructive import make_valid
from shapely.geometry.base import BaseGeometry
from shapely.ops import orient

Expand Down Expand Up @@ -67,7 +68,7 @@ def merge_polygons(polygons: Sequence[BaseGeometry], buffer_distance: float) ->
# Ref: https://datatracker.ietf.org/doc/html/rfc7946#section-3.1.6
oriented_union_simplified = orient(union_rounded, sign=1.0)

return oriented_union_simplified
return make_valid(oriented_union_simplified)


def generate_capture_area(polygons: Sequence[BaseGeometry], gsd: Decimal) -> dict[str, Any]:
Expand Down
23 changes: 23 additions & 0 deletions scripts/stac/imagery/tests/collection_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,13 @@
from moto import mock_aws
from moto.s3.responses import DEFAULT_REGION_NAME
from pytest_subtests import SubTests
from shapely.predicates import is_valid

from scripts.datetimes import format_rfc_3339_datetime_string
from scripts.files.files_helper import ContentType
from scripts.files.fs import read
from scripts.files.fs_s3 import write
from scripts.stac.imagery.capture_area import merge_polygons
from scripts.stac.imagery.collection import ImageryCollection
from scripts.stac.imagery.item import ImageryItem, STACAsset
from scripts.stac.imagery.metadata_constants import CollectionMetadata
Expand Down Expand Up @@ -342,6 +344,27 @@ def test_capture_area_added(fake_collection_metadata: CollectionMetadata, fake_l
assert collection.stac["assets"]["capture_area"]["file:size"] in (269,) # geos 3.11 - geos 3.12 as yet untested


def test_should_make_valid_capture_area() -> None:
# Given two touching triangles
polygons = [
shapely.geometry.shape(
{
"type": "MultiPolygon",
"coordinates": [[[[0, 0], [0, 1], [1, 1], [0, 0]]]],
}
),
shapely.geometry.shape(
{
"type": "MultiPolygon",
"coordinates": [[[[1, 0], [2, 2], [1, 2], [1, 0]]]],
}
),
]

capture_area = merge_polygons(polygons, 0.1)
assert is_valid(capture_area)


def test_event_name_is_present(fake_collection_metadata: CollectionMetadata, fake_linz_slug: str) -> None:
fake_collection_metadata["event_name"] = "Forest Assessment"
collection = ImageryCollection(fake_collection_metadata, any_epoch_datetime, fake_linz_slug)
Expand Down

0 comments on commit be8a039

Please sign in to comment.