Interpolating water depth for multiple scale mesh.

[Jiangchao3]

Describe the bug Hi @krober10nd and @WPringle , I am curious about how does the function of interp() works in interpolating water depth to each vertice for a multiple scale mesh. The result I got seems odd.

To Reproduce Steps to reproduce the behavior: I made a mesh covering two domains, the inner domain is the Pearl River Estuary, the outer domain is the North South China Sea. this is the mesh I generated:

for the inner domain, I use plot(gdat_PRE,'shp') and plot(gdat_PRE,'dem') to plot, get this two figures:

The dem is a rectangularar area. The box of this inner domain I used is an arc line on the side of the sea.

After generating the integrated mesh, I used the interp() function to interpolate water depth to vertices. m_NSCS_PRE = interp(m_NSCS_PRE,{gdat_NSCS,gdat_PRE},'type','depth','nan','fillinside','mindepth',1);

Well, then I used the script plot(m_NSCS_PRE,' type', 'blog'), see the figure below:

The water depth in the log scale around the connected area seems odd. An obvious rectangular area appeared.

So, my question is that when using the interp() function to interpolate water depth to multiple scales mesh, the water depth used for the inner domain seems the whole rectangular dem domain, rather than the domain inside the arc line on the side of the sea.

Expected behavior A clear and concise description of what you expected to happen.

So, I think the water depth interpolated into the inner domain should be bounded by the arc line rather than the whole dem of the inner domain.

I am not sure if my consideration is correct. Please check it.

Screenshots/Figures If applicable, add screenshots to help explain your problem.

Configuration (please complete the following information):

• OS: [e.g. MacOS, Windows, Debian, etc.]
• Version [e.g. OS Mojave]

Additional context Add any other context about the problem here.

[krober10nd]

Hey Jiang

Please don't use interpolation from geodata objects as they downsample the original data for the mesh size function creation. Use the original DEM in NetCDF format.

Secondly, find the indices in your high resolution box and pass them to interp

% bounds of high resolution area 
lon_min  =?? 
lat_min = ?? 
lat_max = ?? 
lon_max = ??
m_NSCS_PRE = interp(m_NSCS_PRE,dem1,'type','depth','nan','fillinside','mindepth',1);
k = find(m_NSCS_PRE.p(:,1) > lon_min &&  m_NSCS_PRE.p(:,1) < lon_max && m_NSCS_PRE.p(:,2) > lat_min && m_NSCS_PRE.p(:,2) < lat_max);   
m_NSCS_PRE = interp(m_NSCS_PRE,dem2,'type','depth','nan','fillinside','mindepth',1,'K',k);

[Jiangchao3]

Hi @krober10nd, Thanks for your reply. However, using the original DEM in NetCDF format file for interpolation just for the outer domain failed because out of the memory.

How should I avoid interpolating out of memory using the original dem NC file?

[krober10nd]

Yes I see. I would recommend looping over the DEM in longitudinal orientated slices using the K index option. Essentially doing the interpolation in sections.

[Jiangchao3]

Well, I will try this looping method you recommended. Seems using the DEM file to interpolation takes a longer run time than using the geodata.

Secondly, I also want to determine one more aspect: when interpolating for the inner high-resolution domain, whether the bounds of the interpolation is the bounds of the original DEM file‘ rectangle ’(because the DEM file is generated in a rectangular bound) rather than the custom arc bound designed by myself?

Because the NC dem file is generated by using a kriging interpolation method. Although the final dem file covers the whole rectangular area, the interpolated scatter data is only within the arc line area. So the water depth value out of the arc line area will be very inaccurate. Look at the followed figure:

Maybe I can use the K index method by setting several lat_min values, these lat_min points can form such an arc line? And any other good suggestion?

Many thanks again!

[Jiangchao3]

Hi @krober10nd , this is my solution:

1、I use the extrac_subdomain function to extract the inner domain and the outer domain, m_NSCS = extract_subdomain(m_NSCS_PRE,bbox_PRE,"keep_inverse",1); plot(m_NSCS,'type','tri') m_PRE = extract_subdomain(m_NSCS_PRE,bbox_PRE); plot(m_PRE,'type','tri');

2、interepolate water depth for the inner domain and the outer domain respectively for j = 1:length(k_index_NSCS) k_lon_min = k_index_NSCS(j,1); k_lon_max = k_index_NSCS(j,2); k_lat_min = k_index_NSCS(j,3); k_lat_max = k_index_NSCS(j,4); k = find(m_NSCS.p(:,1) >= k_lon_min & m_NSCS.p(:,1) < k_lon_max & m_NSCS.p(:,2) > k_lat_min & m_NSCS.p(:,2) < k_lat_max); m_NSCS = interp(m_NSCS,dem_NSCS,'type','depth','nan','fillinside','mindepth',1,'K',k); m_NSCS = interp(m_NSCS,dem_NSCS,'type','slope','slope_calc','rms','K',k); end plot(m_NSCS,'type','blog')

for j = 1:length(k_index_PRE) k_lon_min = k_index_PRE(j,1); k_lon_max = k_index_PRE(j,2); k_lat_min = k_index_PRE(j,3); k_lat_max = k_index_PRE(j,4); k = find(m_PRE.p(:,1) >= k_lon_min & m_PRE.p(:,1) < k_lon_max & m_PRE.p(:,2) > k_lat_min & m_PRE.p(:,2) < k_lat_max); m_PRE = interp(m_PRE,dem_PRE,'type','depth','nan','fillinside','mindepth',1,'K',k); m_PRE = interp(m_PRE,dem_PRE,'type','slope','slope_calc','rms','K',k); end plot(m_PRE,'type','blog');

3、use the plus function to merge these two domain m_NSCS_PRE_new = plus(m_PRE,m_NSCS,'match'); m_NSCS_PRE_new = lim_bathy_slope(m_NSCS_PRE_new,0.1,0); plot(m_NSCS_PRE_new,'type','blog');

This is the first time I know to use the K index option and looping method to interpolate water depth to the vertices. That method will be helpful for making a bigger mesh covering the whole IndWPac domain.

Have to admit that OceanMesh2d is so powerful, and you and William do great works.

I will be grateful if you can give some extra suggestions.

[krober10nd]

Glad you figured it out! I see what your problem was and that was a creative solution!