Use triangle to improve subsetting transition

ocsmesh/cli/subset_n_combine.py

@@ -11,11 +11,12 @@
 from jigsawpy.jig_t import jigsaw_jig_t
 import numpy as np
 from pyproj import CRS, Transformer
+from scipy.interpolate import griddata
 from scipy.spatial import cKDTree
 from shapely.geometry import MultiPolygon, Polygon, GeometryCollection
 from shapely.ops import polygonize, unary_union
 
-from ocsmesh import Raster, Geom, Mesh, Hfun, utils
+from ocsmesh import Raster, Geom, Mesh, Hfun, utils, JigsawDriver
 
 logging.basicConfig(
     stream=sys.stdout,
@@ -242,96 +243,101 @@ def _add_overlap_to_polygon(self, mesh, polygon):
 
         return new_polygon, clipped_mesh
 
-    def _calculate_mesh_size_function(self, hires_mesh_clip, lores_mesh_clip, crs):
+    def _calculate_mesh_size_function(
+            self,
+            buffer_domain,
+            hires_mesh_clip,
+            lores_mesh_clip,
+            buffer_crs
+        ):
 
-        # calculate mesh size for clipped bits
-        hfun_hires = Hfun(Mesh(deepcopy(hires_mesh_clip)))
-        hfun_hires.size_from_mesh()
-        hfun_lowres = Hfun(Mesh(deepcopy(lores_mesh_clip)))
-        hfun_lowres.size_from_mesh()
+        assert(buffer_crs == hires_mesh_clip.crs == lores_mesh_clip.crs)
 
-#    _logger.info("Writing hfun clips...")
-#    start = time()
-#    hfun_hires.mesh.write(str(out_dir / "hfun_fine.2dm"), format="2dm", overwrite=True)
-#    hfun_lowres.mesh.write(str(out_dir / "hfun_coarse.2dm"), format="2dm", overwrite=True)
-#    _logger.info(f"Done in {time() - start} sec")
+        # HARDCODED FOR NOW
+        approx_elem_per_width = 3
 
-        jig_hfun = utils.merge_msh_t(
-            hfun_hires.msh_t(), hfun_lowres.msh_t(),
-            out_crs=crs)#jig_geom) ### TODO: CRS MUST BE == GEOM_MSHT CRS
+        utm = utils.estimate_bounds_utm(
+            buffer_domain.bounds, buffer_crs
+        )
 
-        return jig_hfun
+        # calculate mesh size for clipped bits
+        msht_hi = deepcopy(hires_mesh_clip)
+        utils.reproject(msht_hi, utm)
+        hfun_hi = Hfun(Mesh(msht_hi))
+        hfun_hi.size_from_mesh()
+
+        msht_lo = deepcopy(lores_mesh_clip)
+        utils.reproject(msht_lo, utm)
+        hfun_lo = Hfun(Mesh(msht_lo))
+        hfun_lo.size_from_mesh()
+
+        msht_cdt = utils.triangulate_polygon(
+            buffer_domain, None, opts='p'
+        )
+        msht_cdt.crs = buffer_crs
 
+        utils.reproject(msht_cdt, utm)
+        hfun_cdt = Hfun(Mesh(msht_cdt))
+        hfun_cdt.size_from_mesh()
 
-    def _generate_mesh_for_buffer_region(
-            self, buffer_polygon, jig_hfun, crs):
+        hfun_cdt_sz = deepcopy(hfun_cdt.msh_t().value) / approx_elem_per_width
+        msht_cdt.value[:] = hfun_cdt_sz 
+        hfun_rep = Hfun(Mesh(msht_cdt))
 
-        seam = Geom(buffer_polygon, crs=crs)
+        mesh_domain_rep = JigsawDriver(
+            geom=Geom(buffer_domain, crs=buffer_crs),
+            hfun=hfun_rep,
+            initial_mesh=False
+        ).run()
 
-        jig_geom = seam.msh_t()
+        msht_domain_rep = deepcopy(mesh_domain_rep.msh_t)
+        utils.reproject(msht_domain_rep, utm)
 
-        # IMPORTANT: Setting these to -1 results in NON CONFORMAL boundary
-#    jig_geom.vert2['IDtag'][:] = -1
-#    jig_geom.edge2['IDtag'][:] = -1
+        pts_2mesh = np.vstack(
+            (hfun_hi.msh_t().vert2['coord'], hfun_lo.msh_t().vert2['coord'])
+        )
+        val_2mesh = np.vstack(
+            (hfun_hi.msh_t().value, hfun_lo.msh_t().value)
+        )
+        domain_sz_1 = griddata(
+            pts_2mesh, val_2mesh, msht_domain_rep.vert2['coord'], method='linear'
+        )
+        domain_sz_2 = griddata(
+            pts_2mesh, val_2mesh, msht_domain_rep.vert2['coord'], method='nearest'
+        )
+        domain_sz = domain_sz_1.copy()
+        domain_sz[np.isnan(domain_sz_1)] = domain_sz_2[np.isnan(domain_sz_1)]
 
-        jig_init = deepcopy(seam.msh_t())
-        jig_init.vert2['IDtag'][:] = -1
-        jig_init.edge2['IDtag'][:] = -1
+        msht_domain_rep.value[:] = domain_sz
+        hfun_interp = Hfun(Mesh(msht_domain_rep))
 
-        # Calculate length of all edges on geom
-        geom_edges = jig_geom.edge2['index']
-        geom_coords = jig_geom.vert2['coord'][geom_edges, :]
-        geom_edg_lens = np.sqrt(np.sum(
-            np.power(np.abs(np.diff(geom_coords, axis=1)), 2),
-            axis=2)).squeeze()
+        return hfun_interp
 
-        # TODO: Use marche to calculate mesh size in the area between
-        # the two regions?
 
-        _logger.info("Meshing...")
-        start = time()
-        opts = jigsaw_jig_t()
-        opts.hfun_scal = "absolute"
-        opts.hfun_hmin = np.min(geom_edg_lens)
-        opts.hfun_hmax = np.max(geom_edg_lens)
-#    opts.hfun_hmin = np.min(jig_hfun.value.ravel())
-#    opts.hfun_hmax = np.max(jig_hfun.value.ravel())
-        opts.optm_zip_ = False
-        opts.optm_div_ = False
-        opts.mesh_dims = +2
-        opts.mesh_rad2 = 1.05
-#    opts.mesh_rad2 = 2.0
-
-        jig_mesh = jigsaw_msh_t()
-        jig_mesh.mshID = 'euclidean-mesh'
-        jig_mesh.ndims = 2
-        jig_mesh.crs = jig_geom.crs
-
-        jigsawpy.lib.jigsaw(
-            opts, jig_geom, jig_mesh,
-#        hfun=jig_hfun,
-            init=jig_init
-            )
-
-        jig_mesh.value = np.zeros((jig_mesh.vert2.shape[0], 1))
-        self._transform_mesh(jig_mesh, crs)
-
-        # NOTE: Remove out of domain elements (some corner cases!)
-        elems = jig_mesh.tria3['index']
-        coord = jig_mesh.vert2['coord']
-
-        gdf_elems = gpd.GeoDataFrame(
-            geometry=[Polygon(tri) for tri in coord[elems]],
-            crs=jig_mesh.crs
-        )
-        idx = gdf_elems.representative_point().sindex.query(
-            seam.get_multipolygon(), predicate='intersects'
+    def _generate_mesh_for_buffer_region(
+            self, buffer_polygon, hfun_buffer, buffer_crs):
+
+        mesh_domain_interp = JigsawDriver(
+            geom=Geom(buffer_polygon, crs=buffer_crs),
+            hfun=hfun_buffer,
+            initial_mesh=False
+        ).run()
+
+        msht_domain_interp = deepcopy(mesh_domain_interp.msh_t)
+        in_verts_mask = np.ones_like(msht_domain_interp.vert2, dtype=bool)
+        in_verts_mask[
+            np.unique(utils.get_boundary_edges(msht_domain_interp))
+        ] = False
+
+        utils.reproject(msht_domain_interp, buffer_crs)
+        msht_buffer = utils.triangulate_polygon(
+            shape=buffer_polygon,
+            aux_pts=msht_domain_interp.vert2[in_verts_mask]['coord'],
+            opts='p'
         )
-        flag = np.zeros(len(gdf_elems), dtype=bool)
-        flag[idx] = True
-        jig_mesh.tria3 = jig_mesh.tria3[flag]
+        msht_buffer.crs = buffer_crs
 
-        return jig_mesh
+        return msht_buffer
 
 
     def _transform_mesh(self, mesh, out_crs):
@@ -623,9 +629,13 @@ def _main(
 
         poly_seam = poly_seam_8
 
-        # TODO: Get CRS correctly from geom utm
-#    jig_hfun = self._calculate_mesh_size_function(jig_clip_hires, jig_clip_lowres, crs)
-        jig_buffer_mesh = self._generate_mesh_for_buffer_region(poly_seam, None, crs)
+        hfun_buffer = self._calculate_mesh_size_function(
+            poly_seam, jig_clip_hires, jig_clip_lowres, crs
+        )
+        jig_buffer_mesh = self._generate_mesh_for_buffer_region(
+            poly_seam, hfun_buffer, crs
+        )
+        # TODO: UTM TO CRS?
 
         _logger.info("Combining meshes...")
         start = time()

ocsmesh/geom/base.py

@@ -144,42 +144,8 @@ def multipolygon_to_jigsaw_msh_t(
         transformer = Transformer.from_crs(crs, utm_crs, always_xy=True)
         multipolygon = ops.transform(transformer.transform, multipolygon)
 
-    vert2: List[Tuple[Tuple[float, float], int]] = []
-    for polygon in multipolygon.geoms:
-        if np.all(
-                np.asarray(
-                    polygon.exterior.coords).flatten() == float('inf')):
-            raise NotImplementedError("ellispoidal-mesh")
-        for x, y in polygon.exterior.coords[:-1]:
-            vert2.append(((x, y), 0))
-        for interior in polygon.interiors:
-            for x, y in interior.coords[:-1]:
-                vert2.append(((x, y), 0))
-
-    # edge2
-    edge2: List[Tuple[int, int]] = []
-    for polygon in multipolygon.geoms:
-        polygon = [polygon.exterior, *polygon.interiors]
-        for linear_ring in polygon:
-            _edge2 = []
-            for i in range(len(linear_ring.coords)-2):
-                _edge2.append((i, i+1))
-            _edge2.append((_edge2[-1][1], _edge2[0][0]))
-            edge2.extend(
-                [(e0+len(edge2), e1+len(edge2))
-                    for e0, e1 in _edge2])
-    # geom
-    geom = jigsaw_msh_t()
-    geom.ndims = +2
-    geom.mshID = 'euclidean-mesh'
-    # TODO: Consider ellipsoidal case.
-    # geom.mshID = 'euclidean-mesh' if self._ellipsoid is None \
-    #     else 'ellipsoidal-mesh'
-    geom.vert2 = np.asarray(vert2, dtype=jigsaw_msh_t.VERT2_t)
-    geom.edge2 = np.asarray(
-        [((e0, e1), 0) for e0, e1 in edge2],
-        dtype=jigsaw_msh_t.EDGE2_t)
-    geom.crs = crs
+    msht = utils.shape_to_msh_t(multipolygon)
+    msht.crs = crs
     if utm_crs is not None:
-        geom.crs = utm_crs
-    return geom
+        msht.crs = utm_crs
+    return msht