diff --git a/scripts/python_utilities/coupler/SimGrid.py b/scripts/python_utilities/coupler/SimGrid.py index 5b6a2548..68fe56c9 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"] ####################################### @@ -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) diff --git a/scripts/python_utilities/coupler/couplingUtils.py b/scripts/python_utilities/coupler/couplingUtils.py index 31b84f2a..d2c0d886 100644 --- a/scripts/python_utilities/coupler/couplingUtils.py +++ b/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)) 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 }