Merge branch 'main' into feature/interpolate_other_bands

ocsmesh/hfun/collector.py

@@ -1964,8 +1964,15 @@ def _create_big_raster(self, out_path: Union[str, Path]) -> Raster:
         transformer = Transformer.from_crs(
                 'EPSG:4326', utm_crs, always_xy=True)
 
-        box_epsg4326 = box(x0, y0, x1, y1)
-        poly_utm = ops.transform(transformer.transform, Polygon(box_epsg4326))
+        # If it's the full earth, then the spaces are at least 1/2 deg
+        xs = np.linspace(x0, x1, 720)
+        ys = np.linspace(y0, y1, 720)
+        coords = [[x0, y] for y in ys]
+        coords.extend([[x, y1] for x in xs])
+        coords.extend([[x1, y] for y in reversed(ys)])
+        coords.extend([[x, y0] for x in reversed(xs)])
+        poly_epsg4326 = Polygon(np.array(coords))
+        poly_utm = ops.transform(transformer.transform, poly_epsg4326)
         x0, y0, x1, y1 = poly_utm.bounds
 
         worst_res = 0

tests/api/hfun.py

@@ -461,6 +461,34 @@ def test_add_feature_fast(self):
         self.assertTrue(isinstance(hfun_msht, jigsaw_msh_t))
 
 
+    def test_hfun_fast_extent(self):
+        r_path = self.tdir / 'rast_large.tif'
+        rast_xy = np.mgrid[-80:80:0.2, 20:70:0.2]
+        rast_z = np.ones_like(rast_xy[0]) * 100
+
+        ocsmesh.utils.raster_from_numpy(
+            r_path, rast_z, rast_xy, 4326
+        )
+
+        rast = ocsmesh.Raster(r_path)
+        hfun_coll = ocsmesh.Hfun(
+            [rast], hmin=100000, hmax=200000, method='fast'
+        )
+        hfun_msht = hfun_coll.msh_t()
+
+        ocsmesh.utils.reproject(hfun_msht, rast.crs)
+        rast_box = rast.get_bbox()
+        hfun_box = ocsmesh.utils.get_mesh_polygons(hfun_msht)
+
+        # NOTE: It's good enough if it covers most of it (?)
+        self.assertTrue(
+            rast_box.difference(hfun_box).area / rast_box.area < 5e-4
+        )
+        # This fails due to coarse elements!
+#        self.assertTrue(hfun_box.covers(rast_box)) 
+
+
+
 class SizeFromMesh(unittest.TestCase):
 
     def setUp(self):