I/O routine survived till approx. 15500/21000 grid cells; updated plo…

…tting diagnostics

however, there is a gap between the MERIT and REMA datasets, and the lat-lon grid are also different. This latitudinal strip is not important for orographic data, as only a few small islands exist here, but it is important for me to find a sensible way to glue these two datasets together, e.g., by interpolation.

inputs/icon_regional_run.py

@@ -7,8 +7,6 @@
 params.fn_output = "icon_merit_reg"
 utils.transfer_attributes(params, local_paths.paths, prefix="path")
 
-print(True)
-
 ### alaska
 params.lat_extent = [48.0, 64.0, 64.0]
 params.lon_extent = [-148.0, -148.0, -112.0]
@@ -40,4 +38,5 @@
 params.refine = False
 params.verbose = False
 
-params.plot = True
+params.plot = False
+params.plot_output = True

runs/icon_merit_global.py

@@ -106,8 +106,10 @@ def do_cell(c_idx,
         triangles[i, :, 0] = np.array(clon_vertices[i, :])
         triangles[i, :, 1] = np.array(clat_vertices[i, :])
 
-    if params.plot:
-        cart_plot.lat_lon_icon(topo, triangles, ncells=ncells, clon=clon, clat=clat, title=c_idx)
+    if params.plot or params.plot_output:
+
+        output_fn = params.path_output + str(c_idx) + ".png"
+        cart_plot.lat_lon_icon(topo, triangles, ncells=ncells, clon=clon, clat=clat, title=c_idx, fn = output_fn, output_fig = True)
 
 # %%
     tri_idx = 0
@@ -122,6 +124,7 @@ def do_cell(c_idx,
     sa = interface.second_appx(nhi, nhj, params, topo, tri)
 
     dplot = diagnostics.diag_plotter(params, nhi, nhj)
+    dplot.output_dir = params.path_output
 
 
     tri.tri_lon_verts = triangles[:, :, 0]
@@ -244,9 +247,9 @@ def parallel_wrapper(grid, params, reader, writer):
 
     # lazy_results = []
 
-    # for c_idx in range(n_cells)[:20]:
-    # for c_idx in range(n_cells)[180:190]:
-    for c_idx in range(n_cells)[2046:2060]:
+    for c_idx in range(n_cells)[15455:]:
+    # # for c_idx in range(n_cells)[180:190]:
+    # for c_idx in range(n_cells)[2046:2060]:
         pw_run(c_idx)
     #     lazy_result = dask.delayed(pw_run)(c_idx)
     #     lazy_results.append(lazy_result)

src/io.py

@@ -170,6 +170,7 @@ def __init__(self, cell, params, verbose=False, is_parallel=False):
 
             self.merit_cg = params.merit_cg
             self.split_EW = False
+            self.prev_MERIT = False
 
             if not is_parallel:
                 self.get_topo(cell)
@@ -193,8 +194,12 @@ def get_topo(self, cell):
             lat_min_idx = self.__compute_idx(self.lat_verts.min(), "min", "lat")
             lat_max_idx = self.__compute_idx(self.lat_verts.max(), "max", "lat")
 
-            lon_min_idx = self.__compute_idx(min_lon, "min", "lon")
-            lon_max_idx = self.__compute_idx(max_lon, "max", "lon")
+            if not self.split_EW:
+                lon_min_idx = self.__compute_idx(min_lon, "min", "lon")
+                lon_max_idx = self.__compute_idx(max_lon, "max", "lon")
+            else:
+                lon_min_idx = self.__compute_idx(min_lon, "max", "lon")
+                lon_max_idx = self.__compute_idx(max_lon, "min", "lon")
 
             if ( (self.lon_verts.max() - self.lon_verts.min()) > 180.0 ):
                 lon_idx_rng = list(range(lon_max_idx, len(self.fn_lon) - 1 )) + list(range(0,lon_min_idx + 1))
@@ -227,8 +232,8 @@ def __compute_idx(self, vert, typ, direction):
             if typ == "min":
                 if ((vert - fn_int[where_idx]) < 0.0):
                     if direction == "lon":
-                        if not self.split_EW:
-                            where_idx -= 1
+                        # if not self.split_EW:
+                        where_idx -= 1
                     else:
                         where_idx += 1
             elif typ == "max":
@@ -239,7 +244,7 @@ def __compute_idx(self, vert, typ, direction):
                     else:
                         where_idx -= 1
 
-                if where_idx == (len(fn_int) - 1):
+                if (where_idx == (len(fn_int) - 1)) and self.split_EW:
                     where_idx -= 1
 
             where_idx = int(where_idx)
@@ -445,6 +450,18 @@ def __load_topo(self, cell, fns, dirs, lon_cnt, lat_cnt,  lat_idx_rng, lon_idx_r
                     if n_row == 0:
                         cell.lon += lon[lon_low:lon_high].tolist()
 
+                        # current dataset at n_row = 0 is a MERIT dataset
+                        if "MERIT" in fn:
+                            self.prev_merit = True
+
+                    # topographic data is read over MERIT and REMA interface:
+                    if n_row > 0:
+                        if ("REMA" in fn) and (self.prev_merit):
+                            pass
+
+
+
+
                     lon_sz = lon_high - lon_low
                     lat_sz = lat_high - lat_low
 

vis/cart_plot.py

@@ -427,3 +427,4 @@ def lat_lon_icon(
     # -- maximize and save the PNG file
     if output_fig:
         plt.savefig(fn, bbox_inches="tight", dpi=200)
+        plt.close()