Pre-cull mesh with geojson file

Author: alexolinhager

compass/landice/tests/uummannaq_disko/mesh.py

@@ -1,5 +1,7 @@
+import json
 import os
 
+import jigsawpy
 import numpy as np
 import xarray as xr
 from mpas_tools.scrip.from_mpas import scrip_from_mpas
@@ -90,10 +92,19 @@ def run(self):
 
         logger.info('calling build_cell_width')
 
-        cell_width, x1, y1, geom_points, geom_edges, floodMask = \
-            build_cell_width(
-                self, section_name=section_name,
-                gridded_dataset=source_gridded_dataset_2km)
+        if section_gis.get("define_bnds_by_geojson") == 'True':
+            geom_points, geom_edges = set_geojson_geom_points_and_edges(self)
+
+            cell_width, x1, y1, floodMask = \
+                build_cell_width(self,
+                                 section_name=section_name,
+                                 gridded_dataset=source_gridded_dataset_2km,
+                                 calc_geom_bnds=False)
+        else:
+            cell_width, x1, y1, geom_points, geom_edges, floodMask = \
+                build_cell_width(
+                    self, section_name=section_name,
+                    gridded_dataset=source_gridded_dataset_2km)
 
         build_mali_mesh(
             self, cell_width, x1, y1, geom_points, geom_edges,
@@ -156,3 +167,47 @@ def run(self):
         ds["observedThicknessTendencyUncertainty"] = dHdtErr
         # Write the data to disk
         ds.to_netcdf(self.mesh_filename, 'a')
+
+
+def set_geojson_geom_points_and_edges(self):
+    self.add_input_file(filename='UummannaqDisko.geojson',
+                        package='compass.landice.tests.uummannaq_disko',
+                        target='UummannaqDisko.geojson',
+                        database=None)
+
+    with open('UummannaqDisko.geojson', 'r') as meshMarginFile:
+        geojson_data = json.load(meshMarginFile)
+
+    x = []
+    y = []
+    start_edge = []
+    end_edge = []
+    ct = 0
+
+    for feature in geojson_data['features']:
+        geometry = feature['geometry']
+        for coord in geometry['coordinates']:
+            for sub_coord in coord:
+                x.append(sub_coord[0])
+                y.append(sub_coord[1])
+
+                start_edge.append(ct)
+                end_edge.append(ct + 1)
+                ct = ct + 1
+    x = np.array(x)
+    y = np.array(y)
+    start_edge = np.array(start_edge)
+    end_edge = np.array(end_edge)
+
+    start_edge[-1] = ct - 1
+    end_edge[-1] = 0
+
+    geom_points = np.array([((x[i], y[i]), 0)
+                            for i in range(len(y))],
+                           dtype=jigsawpy.jigsaw_msh_t.VERT2_t)
+
+    geom_edges = np.array([((start_edge[i], end_edge[i]), 0)
+                           for i in range(len(start_edge))],
+                          dtype=jigsawpy.jigsaw_msh_t.EDGE2_t)
+
+    return (geom_points, geom_edges)