Handle DEM as scatter points..

[chlnddev]

Handle DEM as scatter points..

[HamishB]

Hi,

FYI my usual workflow in GRASS is to use:

g.region # to set grid bounds and resolution
r.in.xyz # to import x,y,z data into a raster map. it is very fast and can handle datasets that are many 100s of GB in size and 10s of billions of xyz points. good for lidar and multibeam bathymetry
r.surf.nnbathy # to interpolate over the gaps (an addon module *because it uses the Triangle library which is not GPL compatible)
r.mapcalc # to multiply it by -1 it so depth is positive. can also use the 'zscale' option during the r.in.xyz step
r.out.gdal # to export GeoTIFF

r.contour # to make coastline and maximum elevation limiting contour lines
v.out.ogr # to make shapefiles of those two contour lines

help pages:
https://grass.osgeo.org/grass78/manuals/g.region.html
https://grass.osgeo.org/grass78/manuals/r.in.xyz.html
https://grass.osgeo.org/grass78/manuals/addons/r.surf.nnbathy.html

If working with multibeam bathy data see also MB-Systems. It uses GMT internally to do its gridding and interpolation and can also handle truly massive datasets with grace and ease.

regards,
Hamish

[HamishB]

ps - I wrote this for contour lines -> DEM not xyz -> DEM, but it gives some good examples of how various interpolation techniques compare to each other; if it may be of interest.
https://grasswiki.osgeo.org/wiki/Contour_lines_to_DEM

[WPringle]

Hi @huyquangt

There is no formal way to use xyz data with the "interp" function.

It is very easy to interpolate xyz points to the mesh though:

• scatteredInterpolant:
  F = scatteredInterpolant(x,y,z,'method'); % see help for method options
  m.b = F(m.p);
• Inverse distance weighting (with fixed radius 'rad')
  m.b = IDW(x,y,z,m.p(:,1),m.p(:,2),-2,'fr',rad); %set IDW function attached
  IDW.m.txt

Cheers!

[huyquangtranaus]

Hi Williams, thanks for a prompt reply 👍

Yes, I use scatterIterpolant quite often, but here during mesh generation I want to put constraints dt (timestep, CFL max,min) which require to have bathymetry data read by geodata(varargin)/ParseDEM(obj,varargin)/[xvn, yvn, zvn] = getdemvarnames(fname)...

I could convert scatter points into DEM data and write in a netcdf file, but the file would be very large as I need very a high resolution, so it is impractical for my case.

Is there another thought?

By the way, I still plan to play with your stuff, will keep you updated & stay safe!

Cheers
Huy

[krober10nd]

Hey Huy,

If you want to do multiscale mesh generation around a localized area, the mesh size function can be nested with finer resolution boxes. See examples/Example_ECGC.

[huyquangtranaus]

Hi Keith,

Thanks for that 👍 ... yes what I have been doing is similar to this.

As xyz scatter points cannot be handled in geodata at the present time, I end up with writting several bathymetry files with different resolutions.... so it is fine now..

Cheers,
Huy

[krober10nd]

Hey Huy,

Just to clarify why we did not support xyz data from the beginning of the development of this tool is primarily because we do not want to build a mesh to build a mesh.

In the case of the unstructured xyz data, the only way to use this data for mesh generation is to form a TIN (triangulated irregular network) which is often a Delaunay triangulation. There are several downsides to doing this:

1. Storage of the TIN is memory expensive
2. Query time is significantly slower than a structured grid.
3. Algorithmic complexity is more complex with a TIN. For example, the meshgen/private/smooth_outer.m relies on the topology of the structured grid to do the smoothing of nested sizing functions quickly.

and another point that we considered:

4. Most available data are already DEMs placed on structured grids.

In my experience, survey data (xyz) tend to be less reliable.

And you point out, in your area(s) of interest you could create a structured grid and naively interpolate the xyz to the structured grid creating a DEM in a NETCDF file format. This could be done at several length scales to avoid the memory bottleneck that is inherent to structured grids. And finally, this (memory bottleneck w/ structured grids) is exactly why we developed the mutliscale meshing technique that allows nesting of several sizing function and largely disparate horizontal length scales.

[huyquangtranaus]

Awesome, Keith! Thank for a detailed explanation.. so I got several netcdf-DEM files with different resolutions to serve my multiscale meshes.

However, a "tidy" problem has arisen: when I specify two different regions for 2 different mesh scales, I often see a hole in the overlapping region as attached.

have you seen this before?

[krober10nd]

Hey Huy,

yes, I've seen this before but I've never figured out why exactly it happens, to be honest. One thing that I've suspected in the past was that it was related to msh.clean(). If you could repeat the generation of the mesh without applying cleaning (set cleanup to 0 when calling mshgen) that could help us nail down the source of the bug.

[huyquangtranaus]

Hello Keith,

After playing with grid sizes, I found this issue often occurs when the grid size between two regions (scales) are large.. so my approach is to reduce the different in the grid size between two regions... so the unexpected holes could disappear! Do you want me to close this issue? Have a great weekend, mates!

[krober10nd]

Huy,

That's interesting but did you try turning off the cleaning step? I suspect the sharp mesh size gradient between the nests leads to a poor quality element that's removed for some reason while performing the cleaning step.

Could you just test your examples if you haven't already without cleanup enabled and see what happens (even if you have a steep size gradient between nests)? Thanks and have a great weekend as well.

[huyquangtranaus]

Hi Keith,

what do you mean by turning off the cleaning step? sorry I missed your previous reply :))
does it mean not clear all in the script?

Cheers,
Huy

[krober10nd]

Turn off mesh topology cleaning and error checking by setting the name value argument cleanup --> 0.
mshopts = meshgen('ef',fh,'bou',gdat,'cleanup',0);

[huyquangtranaus]

Thank you Keith,

Yes.. using cleanup = 0, the hole also disappeared!

[krober10nd]

And this was with the original setup that produced the hole?

Could you please send the mesh when you did not apply the cleanup? I'd like to run some tests with it.

[huyquangtranaus]

Hi Keith,

Unfortunately, I can't share files to run my case as there are the licensed bathymetry data for our project attached to it in order to run the mesh. However, as I have seen this issue somewhere before, so I will revisit some meshes I have done which can be shared with you. I will email you when it ready.... Also I could also and some files in examples very soon... Stay safe! Cheers,

[krober10nd]

Hey Huy I don’t need to run the mesh nor the bathymetry. I just the points and face table to figure out what’s wrong with the mesh cleaning routine to produce the error you reported. You can give me the mesh with the bathymetry zeroed out.

On Mon, 29 Jun 2020 at 6:48 PM Huy Quang Tran ***@***.***> wrote: Hi Keith, Unfortunately, I can't share files to run my case as there are the licensed bathymetry data for our project attached to it in order to run the mesh. However, as I have seen this issue somewhere before, so I will revisit some meshes I have done which can be shared with you. I will email you when it ready.... Also I could also and some files in examples very soon... Stay safe! Cheers, — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#91 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEOBZ7CPKQXNMA3EXFLO5DDRZED4PANCNFSM4ODHGVIQ> .

-- -------------------- Keith Roberts Post-doctoral researcher University of São Paulo [email protected] <[email protected]>

[Jiangchao3]

Hey Huy,

Just to clarify why we did not support xyz data from the beginning of the development of this tool is primarily because we do not want to build a mesh to build a mesh.

In the case of the unstructured xyz data, the only way to use this data for mesh generation is to form a TIN (triangulated irregular network) which is often a Delaunay triangulation. There are several downsides to doing this:

1. Storage of the TIN is memory expensive
2. Query time is significantly slower than a structured grid.
3. Algorithmic complexity is more complex with a TIN. For example, the meshgen/private/smooth_outer.m relies on the topology of the structured grid to do the smoothing of nested sizing functions quickly.

and another point that we considered:

1. Most available data are already DEMs placed on structured grids.

In my experience, survey data (xyz) tend to be less reliable.

And you point out, in your area(s) of interest you could create a structured grid and naively interpolate the xyz to the structured grid creating a DEM in a NETCDF file format. This could be done at several length scales to avoid the memory bottleneck that is inherent to structured grids. And finally, this (memory bottleneck w/ structured grids) is exactly why we developed the mutliscale meshing technique that allows nesting of several sizing function and largely disparate horizontal length scales.

Hi @@krober10nd and @huyquangt

Thank you very much for the discussion above, and it provided me with a good idea to solve my problem.

I am trying to build my own mesh with the Oceanmesh2d. I have the same problem as Huy, the bathymetry data I have obtained is also in the form of xyz scatter points. According to Roberts' suggestion, I should creat a structured grid and then interpolate the xyz to the grid creating a DEM in a NETCDF file format. However, there are various software or codes to do this work, as the following website shows: http://www.robertschneiders.de/meshgeneration/software.html#automesh2d. I have no idea to choose.

I want to know which tool can quickly and efficiently generate structured grids, so that I can continue to complete the interpolation and finally use it for the generation of unstructured grids with Oceanmesh2d. Can you recommend it for me? I will appreciate it very much.

Best wishes
Jiangchao

[krober10nd]

So the way I've done in the past is to use grass with a weighted interpolation. See the section about sparse interpolation here https://grasswiki.osgeo.org/wiki/Interpolation then I export to a raster in netcdf format.

Maybe @HamishB can help out as he is one of the programmers of modules in GRASS.

[Jiangchao3]

thanks very much, I have heard of this software, but have not used it. Let me have a try.