FIX: Image tiles crossing the dateline need special handling

When image tile requests crossing a dateline happen, there is
a multipolygon that gets passed to the find images function. We
don't necessarily want just the tiles that are within those
geometries, but also the ones that could be between them.

Author: greglucas

lib/cartopy/io/img_tiles.py

@@ -125,7 +125,40 @@ def _find_images(self, target_domain, target_z, start_tile=(0, 0, 0)):
         # Recursively drill down to the images at the target zoom.
         x0, x1, y0, y1 = self._tileextent(start_tile)
         domain = sgeom.box(x0, y0, x1, y1)
-        if domain.intersects(target_domain):
+
+        # Handle MultiPolygon target domains (e.g., from projections crossing
+        # the International Date Line). For better tile coverage, we check
+        # intersection with any part of the MultiPolygon, and also include
+        # tiles that span between disconnected parts to ensure continuity.
+        intersects = False
+        if hasattr(target_domain, 'geoms'):
+            # MultiPolygon case
+            for geom_part in target_domain.geoms:
+                if domain.intersects(geom_part):
+                    intersects = True
+                    break
+
+            # include tiles that span between disconnected geometry parts
+            # to ensure tiles get selected when crossing the dateline
+            if not intersects and len(target_domain.geoms) > 1:
+                bounds = [geom_part.bounds for geom_part in target_domain.geoms]
+                x_ranges = [(b[0], b[2]) for b in bounds]
+                y_ranges = [(b[1], b[3]) for b in bounds]
+
+                # Check if tile is within the overall bounding range
+                overall_x_min = min(x_range[0] for x_range in x_ranges)
+                overall_x_max = max(x_range[1] for x_range in x_ranges)
+                overall_y_min = min(y_range[0] for y_range in y_ranges)
+                overall_y_max = max(y_range[1] for y_range in y_ranges)
+
+                if (x0 <= overall_x_max and x1 >= overall_x_min and
+                    y0 <= overall_y_max and y1 >= overall_y_min):
+                    intersects = True
+        else:
+            # Single Polygon case
+            intersects = domain.intersects(target_domain)
+
+        if intersects:
             if start_tile[2] == target_z:
                 yield start_tile
             else:

lib/cartopy/tests/test_img_tiles.py

@@ -138,6 +138,33 @@ def test_tile_find_images():
             [(7, 4, 4), (7, 5, 4), (8, 4, 4), (8, 5, 4)])
 
 
+def test_multipolygon_find_images():
+    # Test that _find_images handles MultiPolygon geometries correctly,
+    # particularly for orthographic projections that cross the dateline.
+    gt = cimgt.GoogleTiles()
+
+    # Create a MultiPolygon that simulates the dateline crossing scenario
+    # Two separate rectangles representing the split geometry parts
+    geom1 = sgeom.box(-20037508, -8000000, -12000000, 8000000)  # Left part
+    geom2 = sgeom.box(12000000, -8000000, 20037508, 8000000)   # Right part
+    multi_geom = sgeom.MultiPolygon([geom1, geom2])
+
+    tiles = list(gt.find_images(multi_geom, 2))
+    tile_x_coords = {tile[0] for tile in tiles}
+
+    # Should include all x coordinates from 0 to 3 for complete coverage
+    expected_x_coords = {0, 1, 2, 3}
+    assert tile_x_coords == expected_x_coords
+
+    # Request tiles for the individual pieces
+    tiles_part1 = list(gt.find_images(geom1, 2))
+    tiles_part2 = list(gt.find_images(geom2, 2))
+
+    # Combined individual parts should have fewer tiles than the MultiPolygon
+    # because the MultiPolygon includes intermediate tiles for continuity
+    assert len(tiles) > len(tiles_part1) + len(tiles_part2)
+
+
 @pytest.mark.network
 def test_image_for_domain():
     gt = cimgt.GoogleTiles()