From 372f43db705134310132cac4c76ca25ac5d979e0 Mon Sep 17 00:00:00 2001 From: raluy <126766444+raluy@users.noreply.github.com> Date: Tue, 14 Apr 2026 10:29:27 +0200 Subject: [PATCH 1/6] Block-averaged topography --- scripts/python_utilities/coupler/simgrid.json | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/python_utilities/coupler/simgrid.json b/scripts/python_utilities/coupler/simgrid.json index 02ceee05..5719a16a 100644 --- a/scripts/python_utilities/coupler/simgrid.json +++ b/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 } From 5340a58f695a976635a3d4aeb026f0379b27f161 Mon Sep 17 00:00:00 2001 From: raluy <126766444+raluy@users.noreply.github.com> Date: Tue, 14 Apr 2026 10:31:31 +0200 Subject: [PATCH 2/6] Block-averaged topography --- .../python_utilities/coupler/couplingUtils.py | 34 +++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/scripts/python_utilities/coupler/couplingUtils.py b/scripts/python_utilities/coupler/couplingUtils.py index 31b84f2a..0eeb6afd 100644 --- a/scripts/python_utilities/coupler/couplingUtils.py +++ b/scripts/python_utilities/coupler/couplingUtils.py @@ -45,6 +45,40 @@ 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): + for i in range(Nx): + x0 = x_edges_dst[i] + x1 = x_edges_dst[i + 1] + y0 = y_edges_dst[j] + y1 = y_edges_dst[j + 1] + val = 0.0 + wsum = 0.0 + ix0 = np.searchsorted(x_edges_src, x0, side="right") - 1 + ix1 = np.searchsorted(x_edges_src, x1, side="left") + iy0 = np.searchsorted(y_edges_src, y0, side="right") - 1 + iy1 = np.searchsorted(y_edges_src, y1, side="left") + for iy in range(max(iy0, 0), min(iy1 + 1, ny_src)): + for ix in range(max(ix0, 0), min(ix1 + 1, nx_src)): + xs0 = x_edges_src[ix] + xs1 = x_edges_src[ix + 1] + ys0 = y_edges_src[iy] + ys1 = y_edges_src[iy + 1] + overlap_x = min(x1, xs1) - max(x0, xs0) + overlap_y = min(y1, ys1) - max(y0, ys0) + if overlap_x > 0.0 and overlap_y > 0.0: + w = overlap_x * overlap_y + val += topo_dom_ori[iy, ix] * w + wsum += w + data_topo0[j, i] = val / wsum + return data_topo0 + def smoothTerrain(tPos0,dx): start = time.time() slopeThresh=math.tan(math.radians(35.0)) From 95814dab88cc20c626bdeb5aaa7446ad0a4aa4ee Mon Sep 17 00:00:00 2001 From: raluy <126766444+raluy@users.noreply.github.com> Date: Tue, 14 Apr 2026 10:36:28 +0200 Subject: [PATCH 3/6] Block-averaged topography --- scripts/python_utilities/coupler/SimGrid.py | 26 ++++++++++++++------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/scripts/python_utilities/coupler/SimGrid.py b/scripts/python_utilities/coupler/SimGrid.py index 5b6a2548..15bacb5b 100644 --- a/scripts/python_utilities/coupler/SimGrid.py +++ b/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"] ####################################### @@ -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)) From d25c2591952713473040dee8f5a390182bd1fddf Mon Sep 17 00:00:00 2001 From: raluy <126766444+raluy@users.noreply.github.com> Date: Wed, 15 Apr 2026 10:31:58 +0200 Subject: [PATCH 4/6] Block-averaged topography --- scripts/python_utilities/coupler/SimGrid.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/scripts/python_utilities/coupler/SimGrid.py b/scripts/python_utilities/coupler/SimGrid.py index 15bacb5b..ee90756d 100644 --- a/scripts/python_utilities/coupler/SimGrid.py +++ b/scripts/python_utilities/coupler/SimGrid.py @@ -89,16 +89,14 @@ ## -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)) +npx_inc = d_xi/dx_inter +npy_inc = d_eta/dy_inter +if ((d_xi % dx_inter) == 0.0 and d_xi >= dx_inter): + x_e = x_s + Nx*int(npx_inc) + y_e = y_s + Ny*int(npy_inc) else: - x_e = x_s + Nx*npx_inc - y_e = y_s + Ny*npy_inc + 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] @@ -246,7 +244,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) From c017bda2b2c5c061224f4337c99cad634b1692ec Mon Sep 17 00:00:00 2001 From: raluy <126766444+raluy@users.noreply.github.com> Date: Wed, 15 Apr 2026 10:40:18 +0200 Subject: [PATCH 5/6] Block-averaged topography --- .../python_utilities/coupler/couplingUtils.py | 20 ++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/scripts/python_utilities/coupler/couplingUtils.py b/scripts/python_utilities/coupler/couplingUtils.py index 0eeb6afd..d2c0d886 100644 --- a/scripts/python_utilities/coupler/couplingUtils.py +++ b/scripts/python_utilities/coupler/couplingUtils.py @@ -53,29 +53,31 @@ def block_average_topo(topo_dom_ori,dx_inter,dy_inter,d_xi,d_eta,Nx,Ny): 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] - y0 = y_edges_dst[j] - y1 = y_edges_dst[j + 1] - val = 0.0 - wsum = 0.0 ix0 = np.searchsorted(x_edges_src, x0, side="right") - 1 ix1 = np.searchsorted(x_edges_src, x1, side="left") - iy0 = np.searchsorted(y_edges_src, y0, side="right") - 1 - iy1 = np.searchsorted(y_edges_src, y1, 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] - ys0 = y_edges_src[iy] - ys1 = y_edges_src[iy + 1] overlap_x = min(x1, xs1) - max(x0, xs0) - overlap_y = min(y1, ys1) - max(y0, ys0) 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 From 2a1099fc0354097960ce25f5a95e0c5169e29b0a Mon Sep 17 00:00:00 2001 From: raluy <126766444+raluy@users.noreply.github.com> Date: Wed, 15 Apr 2026 10:58:44 +0200 Subject: [PATCH 6/6] Block-averaged topography --- scripts/python_utilities/coupler/SimGrid.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/scripts/python_utilities/coupler/SimGrid.py b/scripts/python_utilities/coupler/SimGrid.py index ee90756d..68fe56c9 100644 --- a/scripts/python_utilities/coupler/SimGrid.py +++ b/scripts/python_utilities/coupler/SimGrid.py @@ -89,12 +89,14 @@ ## -npx_inc = d_xi/dx_inter -npy_inc = d_eta/dy_inter if ((d_xi % dx_inter) == 0.0 and d_xi >= dx_inter): - x_e = x_s + Nx*int(npx_inc) - y_e = y_s + Ny*int(npy_inc) + 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)