Meshing artefacts at coastline

[acroucher]

I'm creating a mesh with geometry defined by a coastline read in from a shapefile, and the hfun (mesh size function) based on a DEM raster. For the hfun I'm setting the hmin vale to 50 m which is the target resolution at the coast. As in the documentation examples I'm also setting an hfun contour at the high water level with target_size = hmin.

I'm getting some odd mesh artefacts at certain locations along the coast - the resolution is about 50 m along most of it, but in some places it is larger, and in between these there's an area of 50 m elements offshore. The DEM is quite smooth and doesn't seem to have any features that might cause this. Here is a sample plot of the coastline and DEM:

And here is a plot of the hfun:

This gives the following mesh:

Any idea what is causing this to happen?

[SorooshMani-NOAA]

@acroucher this is strange! My first guess is that the DEM might have a concave up bathymetry there (?) which might show up if you change the colormap to only show elevation range smaller than that of the complete figure.
In any case, would it be possible for you to send me the DEM for only the segment shown in the figures above? To do so you can use Raster object and call clip() with the shapely.geometry.bbox of the figures above and then save() the resulting raster. Then you can zip the resulting raster and attach to a comment right here!

[SorooshMani-NOAA]

Looks like there's been no response yet from the Jigsaw dev?

No there wasn't; I think he's very busy with other projects he has on his hands. He was usually very responsive before.

[acroucher]

I'm wondering if another approach to this problem might work better. I tried re-meshing with the hmin reduced from 50 to 25 m, and the larger small islands which had been omitted before were now meshed properly. So it seems Jigsaw is probably deciding these islands are just too small for the 50 m resolution.

Even at 25 m, of course there are yet smaller islands that are omitted, which is OK. I think the problem is really just that there are larger small islands (a few hundred metres across) which I want to be included, and I expected them to be meshed at 50 m but they are not.

However, using 25 m everywhere results in a mesh that is too big and has unnecessary detail in most places. What if I could specify a smaller hmin only for selected islands (which are currently omitted but I want included), and leave hmin = 50 m everywhere else?

Then these islands would be meshed nicely (rather than just clipping out elements), and it would avoid the problems I've been having with other elements being unexpectedly clipped out.

I have had a go at figuring out how to do this, but so far without success - I'm not familiar enough with GeoPandas, OCSMesh etc to get it to work.

At present hmin is applied using a single add_feature():

coast_feature = gpd.read_file(coast_path).to_crs(crs_4326).boundary
hfun.add_feature(coast_feature.unary_union,
                 expansion_rate = 0.004, target_size = hmin)

so I guess this would have to be replaced with some kind of loop over the individual island shapes in coast_feature and doing a separate add_feature() for each one, with the target_size set according to whether the default hmin should be set or if it's one of the special islands where I want higher resolution? Ideally it would be good to be able to identify these islands by name, so I could just have a list of island names (I can see them in QGIS) for the special treatment.

[SorooshMani-NOAA]

I really appreciate your determination in getting this to work with OCSMesh! It's interesting that the smaller element size resolved the issue, I'm glad that you tested it!

That's a very neat idea; I think one way to achieve it would be something like this: The domain Polygon has .interiors property. This refers to the interior holes, i.e. islands. The idea would be to first create the domain geometry, and store its polygon domain_poly = geom.get_multipolygon() then use this polygon to get the island rings:

island_rings = domain_poly.interiors
islands = [Polygon(ring) for ring in island_rings]

or if domain_poly is a MultiPolygon then

island_rings = [isle_ring for poly in domain_poly.geoms for isle_ring in poly.interiors]
islands = [Polygon(ring) for ring in island_rings]

Then find the islands/island rings of interest based either on a threshold for the perimeter (.length) or the area (.area) of the island polygons. Now you can .buffer the resulting filtered islands to get regions where you'd like to have smaller elements.

For the size function, set the minimum to the new 25$m$ that you suggested and then apply the refinements like before (with 50$m$ target size). Now you need to apply a patch refinement for the HfunCollector size function with 25$m$ target size for the regions around the islands:

hfun.add_patch(
  shape=buffered_filtered_islands, # must be in EPGS:4326
  expansion_rate=0.004, # or maybe even `None`
  target_size=25,
 )

And hopefully that's it!

Note that geopandas has similar functionality if you prefer to use it, e.g .interior or .buffer, etc. exists for GeoSeries objects as far as I remember.

[acroucher]

Thanks! this seems to work quite well. Here is the script for it I have tried:

amin, amax = 20e3, 50e3
hmin_refine = 25
buffer_size = 50
domain_poly = geom.get_multipolygon()
island_rings = [isle_ring for poly in domain_poly.geoms for isle_ring in poly.interiors]
islands = [Polygon(ring) for ring in island_rings]
refine_zones = [island.buffer(buffer_size) for island in islands
                        if amin < island.area < amax]
project = pyproj.Transformer.from_crs(geom.crs, crs_4326, always_xy = True)
refine_zones_4326 = MultiPolygon([transform(project.transform, refine_zone)
                                  for refine_zone in refine_zones])
hfun.add_patch(
    shape = refine_zones_4326,
    expansion_rate = 0.006,
    target_size = hmin_refine)

As you can see I have set a minimum island size as well as a maximum one, so I don't apply the refinement to all the little rocks that I don't really want.

[SorooshMani-NOAA]

That's great! Thanks for testing it. Please let me know if you final mesh is acceptable for your use-case or if there's any feature you'd like OCSMesh to have to make life easier!

[acroucher]

I'm also curious to know what happens if you cut DEMs using ocsmesh tools? Can you please try that as well?

I just tried this, but it doesn't seem to be getting the hfun quite right on the eastern side of the open sea boundary. There is unnecessary refinement there (this also happens if I generate the hfun with the hi-res DEM only). There are also some glitches along the boundary between the hi-res and low-res DEMs. Plus, trying to interpolate the bathymetry from the hi-res and clipped low-res rasters this way does not work at all - I get mesh elevations set to -99999 outside of the hi-res DEM.

[SorooshMani-NOAA]

OK, thanks for checking, I will reword the documentation (#40). If you get a chance to test it with inverse shape, please let me know if it resolves the issue and obviates the need for external tools for cutting the DEMs.

[acroucher]

Is there an easy way to create the inverse shape? (Currently the shape is the bbox of the hi-res DEM raster.)

[SorooshMani-NOAA]

The easiest way I can think of is to get bbox of the low-res DEM and subtract hires bbox from it to get the inverse shape. You can call get_bbox() to get the bounding box of a raster object:

OCSMesh/ocsmesh/raster.py

Lines 605 to 609 in a105e9a

	def get_bbox(
	self,
	crs: Union[str, CRS, None] = None,
	output_type: Literal['polygon', 'bbox'] = 'polygon'
	) -> Union[Polygon, Bbox]:

This should give you a Shapely Polygon object that you can use to calculate the inverse shape.

[acroucher]

OK, I have tried this. It does work, but the hi-res and clipped low-res DEMs don't necessarily quite overlap - I guess perhaps because they have different resolutions. So there could be a gap between the two DEMs, which I figured was probably not good - might lead to glitches along the boundary. So I buffered the clip shape by 100 m to make sure they overlap slightly:

nz_dem = Raster(nz_dem_path)
hauraki_dem = Raster(hauraki_dem_path)
clip_shape = nz_dem.get_bbox().difference(hauraki_dem.get_bbox())
clip_shape = clip_shape.buffer(100, single_sided = True)
nz_dem.clip(clip_shape)

This works and produces a mesh with no glitches and has the correct resolution on the open sea boundary.

[SorooshMani-NOAA]

I see, I've also had my share of issues with bbox overlap for DEMs. Sometimes CRS plays a role too. But the good news is with this approach (with the added buffer) you'll have an automated way all the way from DEM+shapes to mesh without the need for external tools. If you see any opportunity for improvement in the bbox functionality, please let me know. Thanks!