Skip to content

Commit 2b5a2a6

Browse files
committed
fix(affine): protect use_affine and support convert
1 parent acc053e commit 2b5a2a6

File tree

2 files changed

+365
-4
lines changed

2 files changed

+365
-4
lines changed

src/easyidp/geotiff.py

Lines changed: 191 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -299,6 +299,194 @@ def set_mask_polygon(
299299
self._mask_polygon_is_geo = is_geo
300300
# Clear cached binary mask when polygon changes
301301
self._mask = None
302+
303+
def convert_to_affine(self) -> 'GeoTiff':
304+
"""Convert from standard storage to affine rotation storage.
305+
306+
Returns a NEW GeoTiff object with imarray aligned to the mask polygon
307+
rectangle, and transform including rotation. The original object is
308+
not modified.
309+
310+
Requires mask_polygon to be a valid rectangle (4 vertices, 90° angles).
311+
312+
Returns
313+
-------
314+
GeoTiff
315+
New GeoTiff object in affine mode. Returns self if already affine.
316+
317+
Raises
318+
------
319+
ValueError
320+
If no mask polygon is set or polygon is not a valid rectangle.
321+
322+
Example
323+
-------
324+
.. code-block:: python
325+
326+
>>> gtiff = idp.GeoTiff('input.tif')
327+
>>> gtiff.set_mask_polygon(rect_coords, is_geo=True)
328+
>>> affine_gtiff = gtiff.convert_to_affine()
329+
>>> affine_gtiff.use_affine
330+
True
331+
>>> gtiff.use_affine # Original unchanged
332+
False
333+
"""
334+
if self._use_affine:
335+
logger.warning("Already in affine mode, returning self")
336+
return self
337+
338+
if self._mask_polygon is None:
339+
raise ValueError("No mask polygon set. Use set_mask_polygon first.")
340+
341+
polygon_geo = self.mask_polygon_geo
342+
is_rect, angle, bounds = self._is_valid_rectangle(polygon_geo)
343+
344+
if not is_rect:
345+
raise ValueError(
346+
"Polygon is not a valid rectangle. Cannot convert to affine mode."
347+
)
348+
349+
# Perform the transformation
350+
profile = self.header['profile'].copy()
351+
new_imarray, new_profile = self._prepare_affine_storage(
352+
self._imarray.copy(), profile, polygon_geo, angle, bounds
353+
)
354+
355+
# Build new header
356+
new_header = self.header.copy()
357+
new_header['profile'] = new_profile
358+
new_header['width'] = new_profile['width']
359+
new_header['height'] = new_profile['height']
360+
new_header['transform'] = new_profile['transform']
361+
new_header['scale'] = [abs(new_profile['transform'].a),
362+
abs(new_profile['transform'].e)]
363+
364+
# Create new GeoTiff object
365+
new_gtiff = GeoTiff(imarray=new_imarray, header=new_header)
366+
new_gtiff._mask_polygon = self._mask_polygon.copy()
367+
new_gtiff._mask_polygon_is_geo = self._mask_polygon_is_geo
368+
new_gtiff._use_affine = True
369+
370+
logger.info(f"Converted to affine mode with {angle:.1f}° rotation")
371+
return new_gtiff
372+
373+
def convert_from_affine(self, target_bounds: tuple | None = None) -> 'GeoTiff':
374+
"""Convert from affine rotation storage to standard storage.
375+
376+
Returns a NEW GeoTiff object with the rotated imarray transformed back
377+
to axis-aligned coordinates, with a mask for the valid region.
378+
The original object is not modified.
379+
380+
Parameters
381+
----------
382+
target_bounds : tuple, optional
383+
(min_x, min_y, max_x, max_y) for the output image bounds.
384+
If None, uses the bounding box of mask_polygon.
385+
386+
Returns
387+
-------
388+
GeoTiff
389+
New GeoTiff object in standard mode. Returns self if already standard.
390+
391+
Raises
392+
------
393+
ValueError
394+
If no mask polygon is available for conversion.
395+
396+
Example
397+
-------
398+
.. code-block:: python
399+
400+
>>> gtiff = idp.GeoTiff('affine_file.tif')
401+
>>> gtiff.use_affine
402+
True
403+
>>> standard_gtiff = gtiff.convert_from_affine()
404+
>>> standard_gtiff.use_affine
405+
False
406+
>>> gtiff.use_affine # Original unchanged
407+
True
408+
"""
409+
if not self._use_affine:
410+
logger.warning("Already in standard mode, returning self")
411+
return self
412+
413+
if self._mask_polygon is None:
414+
raise ValueError("No mask polygon available for conversion.")
415+
416+
from rasterio.transform import Affine
417+
from scipy.ndimage import map_coordinates
418+
419+
polygon_geo = self.mask_polygon_geo
420+
old_transform = self.header['profile']['transform']
421+
422+
# Calculate target bounds
423+
if target_bounds is None:
424+
min_x, max_x = polygon_geo[:, 0].min(), polygon_geo[:, 0].max()
425+
min_y, max_y = polygon_geo[:, 1].min(), polygon_geo[:, 1].max()
426+
else:
427+
min_x, min_y, max_x, max_y = target_bounds
428+
429+
# Use same pixel size
430+
pixel_size_x = abs(old_transform.a)
431+
pixel_size_y = abs(old_transform.e)
432+
433+
# Calculate output dimensions
434+
out_width = int(np.ceil((max_x - min_x) / pixel_size_x))
435+
out_height = int(np.ceil((max_y - min_y) / pixel_size_y))
436+
437+
# New axis-aligned transform
438+
new_transform = Affine.translation(min_x, max_y) * Affine.scale(pixel_size_x, -pixel_size_y)
439+
440+
# Create output coordinate grids
441+
out_rows, out_cols = np.mgrid[0:out_height, 0:out_width]
442+
443+
# Output pixels to geo coordinates (axis-aligned)
444+
geo_x = min_x + out_cols * pixel_size_x
445+
geo_y = max_y - out_rows * pixel_size_y
446+
447+
# Geo coordinates to input (rotated) pixel coordinates
448+
inv_old = ~old_transform
449+
in_cols = inv_old.a * geo_x + inv_old.b * geo_y + inv_old.c
450+
in_rows = inv_old.d * geo_x + inv_old.e * geo_y + inv_old.f
451+
452+
# Sample from rotated image
453+
imarray = self._imarray
454+
ndim = len(imarray.shape)
455+
if ndim == 2:
456+
new_imarray = map_coordinates(
457+
imarray, [in_rows, in_cols], order=1, mode='constant', cval=0
458+
).astype(imarray.dtype)
459+
else:
460+
bands = imarray.shape[2]
461+
new_imarray = np.zeros((out_height, out_width, bands), dtype=imarray.dtype)
462+
for b in range(bands):
463+
new_imarray[:, :, b] = map_coordinates(
464+
imarray[:, :, b], [in_rows, in_cols],
465+
order=1, mode='constant', cval=0
466+
).astype(imarray.dtype)
467+
468+
# Build new header
469+
new_profile = self.header['profile'].copy()
470+
new_profile['width'] = out_width
471+
new_profile['height'] = out_height
472+
new_profile['transform'] = new_transform
473+
474+
new_header = self.header.copy()
475+
new_header['profile'] = new_profile
476+
new_header['width'] = out_width
477+
new_header['height'] = out_height
478+
new_header['transform'] = new_transform
479+
new_header['scale'] = [pixel_size_x, pixel_size_y]
480+
new_header['tie_point'] = [min_x, max_y]
481+
482+
# Create new GeoTiff object
483+
new_gtiff = GeoTiff(imarray=new_imarray, header=new_header)
484+
new_gtiff._mask_polygon = self._mask_polygon.copy()
485+
new_gtiff._mask_polygon_is_geo = self._mask_polygon_is_geo
486+
new_gtiff._use_affine = False
487+
488+
logger.info("Converted from affine to standard mode")
489+
return new_gtiff
302490

303491
def _polygon_pixel_to_geo(self, polygon: np.ndarray) -> np.ndarray:
304492
"""Convert pixel polygon to geo coordinates using transform."""
@@ -743,12 +931,11 @@ def save(
743931
if use_affine:
744932
is_rect, angle, bounds = self._is_valid_rectangle(polygon_geo)
745933
if is_rect:
746-
# Prepare affine rotated storage
934+
# Prepare affine rotated storage (only for saved file)
747935
imarray, profile = self._prepare_affine_storage(
748936
imarray, profile, polygon_geo, angle, bounds
749937
)
750938
logger.info(f"Affine mode: saved with {angle:.1f}° rotation")
751-
self._use_affine = True
752939
apply_mask = False # No mask needed in affine mode
753940
else:
754941
logger.warning(
@@ -811,8 +998,8 @@ def _is_valid_rectangle(
811998
bounds : tuple
812999
(origin_x, origin_y, width, height) of the rectangle.
8131000
"""
814-
# Remove closure point if present
815-
pts = polygon[:-1] if np.allclose(polygon[0], polygon[-1]) else polygon
1001+
# Remove closure point if present (use strict atol for geo coords)
1002+
pts = polygon[:-1] if np.allclose(polygon[0], polygon[-1], rtol=0, atol=1e-6) else polygon
8161003

8171004
if len(pts) != 4:
8181005
return False, 0.0, ()

tests/test_geotiff.py

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1058,3 +1058,177 @@ def test_one_raw_roi2geotiff_stores_polygon(self, shared_data):
10581058
assert gtiff.mask_polygon is not None
10591059
# Polygon should match input (approximately closed)
10601060
assert len(gtiff.mask_polygon) >= len(roi_geo) - 1
1061+
1062+
1063+
class TestAffineConversion:
1064+
"""Tests for affine mode conversion functions."""
1065+
1066+
def test_convert_to_affine_returns_new_object(self, shared_data):
1067+
"""convert_to_affine returns a new GeoTiff, original unchanged."""
1068+
test_data = shared_data['test_data']
1069+
1070+
# Load and set polygon
1071+
gtiff = idp.GeoTiff(test_data.pix4d.lotus_dom_part)
1072+
gtiff._imarray = gtiff.imarray # Force load
1073+
original_shape = gtiff.imarray.shape
1074+
1075+
# Set axis-aligned rectangle
1076+
polygon = np.array([
1077+
[368025.0, 3955479.0],
1078+
[368027.0, 3955479.0],
1079+
[368027.0, 3955477.0],
1080+
[368025.0, 3955477.0],
1081+
])
1082+
gtiff.set_mask_polygon(polygon, is_geo=True)
1083+
1084+
# Convert to affine
1085+
affine_gtiff = gtiff.convert_to_affine()
1086+
1087+
# Verify new object
1088+
assert affine_gtiff is not gtiff
1089+
assert affine_gtiff.use_affine is True
1090+
1091+
# Verify original unchanged
1092+
assert gtiff.use_affine is False
1093+
assert gtiff.imarray.shape == original_shape
1094+
1095+
def test_convert_from_affine_returns_new_object(self, shared_data):
1096+
"""convert_from_affine returns a new GeoTiff, original unchanged."""
1097+
test_data = shared_data['test_data']
1098+
1099+
# Load, set polygon, convert to affine
1100+
gtiff = idp.GeoTiff(test_data.pix4d.lotus_dom_part)
1101+
gtiff._imarray = gtiff.imarray
1102+
1103+
polygon = np.array([
1104+
[368025.0, 3955479.0],
1105+
[368027.0, 3955479.0],
1106+
[368027.0, 3955477.0],
1107+
[368025.0, 3955477.0],
1108+
])
1109+
gtiff.set_mask_polygon(polygon, is_geo=True)
1110+
1111+
affine_gtiff = gtiff.convert_to_affine()
1112+
affine_shape = affine_gtiff.imarray.shape
1113+
1114+
# Convert back from affine
1115+
standard_gtiff = affine_gtiff.convert_from_affine()
1116+
1117+
# Verify new object
1118+
assert standard_gtiff is not affine_gtiff
1119+
assert standard_gtiff.use_affine is False
1120+
1121+
# Verify affine original unchanged
1122+
assert affine_gtiff.use_affine is True
1123+
assert affine_gtiff.imarray.shape == affine_shape
1124+
1125+
def test_convert_roundtrip_data_consistency(self, shared_data):
1126+
"""Roundtrip conversion preserves polygon and basic properties."""
1127+
test_data = shared_data['test_data']
1128+
1129+
# Load original image
1130+
gtiff = idp.GeoTiff(test_data.pix4d.lotus_dom_part)
1131+
gtiff._imarray = gtiff.imarray
1132+
original_shape = gtiff.imarray.shape
1133+
1134+
# Set axis-aligned rectangle
1135+
polygon = np.array([
1136+
[368025.0, 3955479.0],
1137+
[368027.0, 3955479.0],
1138+
[368027.0, 3955477.0],
1139+
[368025.0, 3955477.0],
1140+
])
1141+
gtiff.set_mask_polygon(polygon, is_geo=True)
1142+
1143+
# Convert to affine
1144+
affine_gtiff = gtiff.convert_to_affine()
1145+
assert affine_gtiff.use_affine is True
1146+
# Affine should be smaller (cropped to rectangle)
1147+
assert affine_gtiff.width <= gtiff.width
1148+
assert affine_gtiff.height <= gtiff.height
1149+
1150+
# Convert back
1151+
roundtrip_gtiff = affine_gtiff.convert_from_affine()
1152+
assert roundtrip_gtiff.use_affine is False
1153+
1154+
# Polygon should be preserved through both conversions
1155+
assert roundtrip_gtiff.mask_polygon is not None
1156+
np.testing.assert_allclose(
1157+
roundtrip_gtiff.mask_polygon_geo[:4],
1158+
polygon,
1159+
atol=1e-6
1160+
)
1161+
1162+
# Roundtrip should have data (not all zeros)
1163+
assert roundtrip_gtiff.imarray.sum() > 0
1164+
1165+
def test_convert_preserves_original(self, shared_data):
1166+
"""Verify original GeoTiff is completely unchanged after conversion."""
1167+
test_data = shared_data['test_data']
1168+
1169+
gtiff = idp.GeoTiff(test_data.pix4d.lotus_dom_part)
1170+
gtiff._imarray = gtiff.imarray
1171+
1172+
# Store original state
1173+
original_imarray = gtiff.imarray.copy()
1174+
original_use_affine = gtiff.use_affine
1175+
original_width = gtiff.width
1176+
original_height = gtiff.height
1177+
1178+
polygon = np.array([
1179+
[368025.0, 3955479.0],
1180+
[368027.0, 3955479.0],
1181+
[368027.0, 3955477.0],
1182+
[368025.0, 3955477.0],
1183+
])
1184+
gtiff.set_mask_polygon(polygon, is_geo=True)
1185+
1186+
# Convert to affine
1187+
_ = gtiff.convert_to_affine()
1188+
1189+
# Verify original is unchanged
1190+
assert gtiff.use_affine == original_use_affine
1191+
assert gtiff.width == original_width
1192+
assert gtiff.height == original_height
1193+
np.testing.assert_array_equal(gtiff.imarray, original_imarray)
1194+
1195+
def test_convert_non_rectangle_raises(self, shared_data):
1196+
"""convert_to_affine raises error for non-rectangle polygons."""
1197+
test_data = shared_data['test_data']
1198+
1199+
gtiff = idp.GeoTiff(test_data.pix4d.lotus_dom_part)
1200+
gtiff._imarray = gtiff.imarray
1201+
1202+
# Set triangle (not rectangle)
1203+
triangle = np.array([
1204+
[368025.0, 3955479.0],
1205+
[368027.0, 3955479.0],
1206+
[368026.0, 3955477.0],
1207+
])
1208+
gtiff.set_mask_polygon(triangle, is_geo=True)
1209+
1210+
# Should raise ValueError
1211+
with pytest.raises(ValueError, match="not a valid rectangle"):
1212+
gtiff.convert_to_affine()
1213+
1214+
def test_already_affine_returns_self(self, shared_data):
1215+
"""convert_to_affine returns self if already in affine mode."""
1216+
test_data = shared_data['test_data']
1217+
1218+
gtiff = idp.GeoTiff(test_data.pix4d.lotus_dom_part)
1219+
gtiff._imarray = gtiff.imarray
1220+
1221+
polygon = np.array([
1222+
[368025.0, 3955479.0],
1223+
[368027.0, 3955479.0],
1224+
[368027.0, 3955477.0],
1225+
[368025.0, 3955477.0],
1226+
])
1227+
gtiff.set_mask_polygon(polygon, is_geo=True)
1228+
1229+
# Convert twice
1230+
affine1 = gtiff.convert_to_affine()
1231+
affine2 = affine1.convert_to_affine()
1232+
1233+
# Second call should return same object
1234+
assert affine2 is affine1

0 commit comments

Comments
 (0)