Eloisa's block-averaged topography and some minor fixes

scripts/python_utilities/coupler/SimGrid.py

@@ -11,6 +11,7 @@
 import math
 import struct
 from scipy.interpolate import RectBivariateSpline, NearestNDInterpolator
+from skimage.measure import block_reduce
 from couplingUtils import *
 
 #######################################
@@ -31,6 +32,7 @@
 name_dom_add = params["name_dom_add"]
 urban_heatRedis_opt = params["urban_heatRedis_opt"]
 landcover_table = params["landcover_table"]
+topo_average_opt = params["topo_average_opt"]
 save_plot_opt = params["save_plot_opt"]
 
 #######################################
@@ -87,16 +89,16 @@
 
 ##
 
-npx_inc = int(d_xi/dx_inter)
-npy_inc = int(d_eta/dy_inter)
-if (npx_inc==0):
-    npx_inc = d_xi/dx_inter
-    npy_inc = d_eta/dy_inter
-    x_e = x_s + int(np.ceil(Nx*npx_inc))
-    y_e = y_s + int(np.ceil(Ny*npy_inc))
-else:
+if ((d_xi % dx_inter) == 0.0 and d_xi >= dx_inter):
+    npx_inc = int(d_xi/dx_inter)
+    npy_inc = int(d_eta/dy_inter)
     x_e = x_s + Nx*npx_inc
     y_e = y_s + Ny*npy_inc
+else:
+    npx_inc = d_xi/dx_inter
+    npy_inc = d_eta/dy_inter
+    x_e = x_s + int(np.ceil((Nx-1)*npx_inc)) + 1
+    y_e = y_s + int(np.ceil((Ny-1)*npy_inc)) + 1
 print('x_s,x_e,y_s,y_e=',x_s,x_e,y_s,y_e)
 
 box_indx = [x_s,x_e,x_e,x_s,x_s]
@@ -115,7 +117,9 @@
     print('Grid is not an even factor of GIS resolution, use interpolation')
     interp_flag = 1
 print('interp_flag=',interp_flag)
-
+if (topo_average_opt==1):
+    print('Using block averaging for topography')
+
 verticalDeformSwitch = int(str(FE_params['verticalDeformSwitch'][0]))
 print('verticalDeformSwitch=',verticalDeformSwitch)
 if (verticalDeformSwitch==1):
@@ -132,24 +136,28 @@
 ## Read in terrain elevation array
 topo = ds_GIS.topoPos.values
 if (interp_flag==0):
-    data_topo0 = topo[y_s:y_e:npy_inc,x_s:x_e:npx_inc]
+    if (topo_average_opt==0):
+        data_topo0 = topo[y_s:y_e:npy_inc,x_s:x_e:npx_inc]
+    else:
+        data_topo0 = block_reduce(topo[y_s:y_e,x_s:x_e], block_size=(int(d_eta/dy_inter),int(d_xi/dx_inter)), func=np.mean)
 else:
     xPos_2d_dom_ori = xPos_2d[y_s:y_e,x_s:x_e]
     yPos_2d_dom_ori = yPos_2d[y_s:y_e,x_s:x_e]
-    topo_dom_ori = topo[y_s:y_e,x_s:x_e]
     print('xPos_2d_dom_ori.shape=',xPos_2d_dom_ori.shape)
-    f_topo = RectBivariateSpline(xPos_2d_dom_ori[0,:], yPos_2d_dom_ori[:,0], topo_dom_ori.T, kx=3, ky=3)
-
     xPos_1d_new = np.arange(x_box_corners[0],x_box_corners[1],d_xi)
     yPos_1d_new = np.arange(y_box_corners[0],y_box_corners[2],d_eta)
-
-    data_topo0_b = f_topo(xPos_1d_new, yPos_1d_new).T
     xPos_2d_new_b, yPos_2d_new_b = np.meshgrid(xPos_1d_new, yPos_1d_new)
-
-    data_topo0 = data_topo0_b[0:Ny,0:Nx]
     xPos_2d_new = xPos_2d_new_b[0:Ny,0:Nx]
     yPos_2d_new = yPos_2d_new_b[0:Ny,0:Nx]
 
+    topo_dom_ori = topo[y_s:y_e,x_s:x_e]
+    if (topo_average_opt==0):
+        f_topo = RectBivariateSpline(xPos_2d_dom_ori[0,:], yPos_2d_dom_ori[:,0], topo_dom_ori.T, kx=3, ky=3)
+        data_topo0_b = f_topo(xPos_1d_new, yPos_1d_new).T
+        data_topo0 = data_topo0_b[0:Ny,0:Nx]
+    else:
+        data_topo0 = block_average_topo(topo_dom_ori,dx_inter,dy_inter,d_xi,d_eta,Nx,Ny)
+
 data_topo = smoothTerrain(data_topo0,d_xi)
 
 topoPos_min = np.min(data_topo,axis=(1,0))
@@ -238,7 +246,7 @@
     if (interp_flag==0):
         data_bmask = bdg_heights[y_s:y_e:npy_inc,x_s:x_e:npx_inc]
     else:
-        f_bdg = NearestNDInterpolator(list(zip(xPos_2d_dom_ori.flatten(), yPos_2d_dom_ori.flatten())), data_bmask[y_s:y_e,x_s:x_e].flatten())
+        f_bdg = NearestNDInterpolator(list(zip(xPos_2d_dom_ori.flatten(), yPos_2d_dom_ori.flatten())), bdg_heights[y_s:y_e,x_s:x_e].flatten())
         data_bmask = f_bdg(xPos_2d_new, yPos_2d_new)
 
     bdg3d_tmp = np.zeros((Nz,Ny,Nx),dtype=np.float32)

scripts/python_utilities/coupler/couplingUtils.py

@@ -45,6 +45,42 @@ def zDeform(zRect, zGround, zCeiling, c1, fCoeff):
 
     return zStretch
 
+def block_average_topo(topo_dom_ori,dx_inter,dy_inter,d_xi,d_eta,Nx,Ny):
+    ny_src, nx_src = topo_dom_ori.shape
+    data_topo0 = np.zeros((Ny, Nx), dtype=float)
+    x_edges_src = np.arange(nx_src + 1) * dx_inter
+    y_edges_src = np.arange(ny_src + 1) * dy_inter
+    x_edges_dst = np.arange(Nx + 1) * d_xi
+    y_edges_dst = np.arange(Ny + 1) * d_eta
+    for j in range(Ny):
+        y0 = y_edges_dst[j]
+        y1 = y_edges_dst[j + 1]
+        iy0 = np.searchsorted(y_edges_src, y0, side="right") - 1
+        iy1 = np.searchsorted(y_edges_src, y1, side="left")
+        for i in range(Nx):
+            x0 = x_edges_dst[i]
+            x1 = x_edges_dst[i + 1]
+            ix0 = np.searchsorted(x_edges_src, x0, side="right") - 1
+            ix1 = np.searchsorted(x_edges_src, x1, side="left")
+            val = 0.0
+            wsum = 0.0
+            for iy in range(max(iy0, 0), min(iy1 + 1, ny_src)):
+                ys0 = y_edges_src[iy]
+                ys1 = y_edges_src[iy + 1]
+                overlap_y = min(y1, ys1) - max(y0, ys0)
+                for ix in range(max(ix0, 0), min(ix1 + 1, nx_src)):
+                    xs0 = x_edges_src[ix]
+                    xs1 = x_edges_src[ix + 1]
+                    overlap_x = min(x1, xs1) - max(x0, xs0)
+                    if overlap_x > 0.0 and overlap_y > 0.0:
+                        w = overlap_x * overlap_y
+                        val += topo_dom_ori[iy, ix] * w
+                        wsum += w
+            if wsum == 0.0:
+                raise ValueError(f'No overlap found for cell j={j}, i={i}.')
+            data_topo0[j, i] = val / wsum
+    return data_topo0
+
 def smoothTerrain(tPos0,dx):
     start = time.time()
     slopeThresh=math.tan(math.radians(35.0))

scripts/python_utilities/coupler/simgrid.json

@@ -9,5 +9,6 @@
   "name_dom_add": "",
   "urban_heatRedis_opt": 0,
   "landcover_table": "/path/to/LandCoverMetadata_NLCD16.csv",
+  "topo_average_opt": 0,
   "save_plot_opt": 1
 }