Masking nodata when reading in a DEM file?

[melissamattwig]

Is there a way to mask nodata values in a DEM? I'm trying to create a large mesh of Lake Ontario and I'm starting with using a very large DEM to start before adding additional, higher resolution lidar data. I added two contours just to see what it would look like (on the left) and also tried plotting the values of the original DEM raster file and I'm getting a very binary looking plot (on the right).

When I made a plot of the histogram for this DEM file it seems like there is a lot of nodata or fill values because when I read in the DEM using rasterio and having masked=True, the values are a lot more evenly distributed.

Please let me know if there is a solution or if I am misinterpreting how I should be working with this DEM. Thank you!

[SorooshMani-NOAA]

As you said the right plot could be just due to large fill values, but I don't understand what the problem is with that plot in terms of OCSMesh use?

For the size function, if your concern is having a lot of elements in the regions where you don't have any data, I suggest using a large hmax, so that all the regions without actual bathymetry data will get the largest size. The raster size function takes the whole DEM bbox for meshing. Then add a contour or constant value size of your actual hmax size at the max elevation level you're interested in.

[melissamattwig]

The plot on the right is of the ocsmesh.raster.Raster object before I add any contours which I was hoping to use as a guide to determine where to add contours based on the bathymetry. I'd never used rasterio before starting to work with ocsmesh, so I'm having a hard time plotting the raster and masking (or at least just not plotting) all of the large fill values that are added within the bounding box of the raster but outside of the actual DEM data.

I was, however, able to look at the DEM in ArcGIS Pro and tried adding the contour lines based on that. It seems like the way that the contours ended up is opposite of what I intended. I wanted the smaller element sizes along the shoreline/floodplain of the lake and then increase towards the center, but it looks like the smaller element sizes are in the middle and there are only a couple contours (I added 8 in total). The code I have for adding the contours is below:

# set value for mean surface level for Lake Ontario
ref_r1 = 74.2 ##coastline

## add contours and constraints
print("adding contours")
print('adding topo bound constraint')
hfun_r1.add_topo_bound_constraint(value=40, upper_bound=ref_r1+50, lower_bound=ref_r1-50, value_type='min', rate=0.01)
print('floodplain')
hfun_r1.add_contour(level=ref_r1+50,expansion_rate=0.0005, target_size=150)
print('shoreline')
hfun_r1.add_contour(level=ref_r1,expansion_rate=0.0005, target_size=75)
print('add constant value between floodplain and shoreline')
hfun_r1.add_constant_value(value=200, lower_bound=ref_r1, upper_bound=ref_r1+50)

hfun_r1.add_contour(level=ref_r1-75,expansion_rate=0.0005, target_size=300)
hfun_r1.add_contour(level=ref_r1-150,expansion_rate=0.0005, target_size=800)
hfun_r1.add_contour(level=ref_r1-200,expansion_rate=0.001, target_size=1500)
hfun_r1.add_contour(level=ref_r1-225,expansion_rate=0.001, target_size=2000)
hfun_r1.add_contour(level=ref_r1-245,expansion_rate=0.001, target_size=2700)

The ref_r1 value is 74.5 which is the mean surface elevation for Lake Ontario and generally where the coastline is and is also what I used to reference all the levels for the contours. In ArcGIS Pro the DEM got more negative deeper into the lake and more positive outside of the lake which is how I approached adding the contours but that's not how it ended up in the size function (pictured below on the left).

It looks a little better than the first attempt I made which I pictured in the original discussion comment I made, but I was trying to get it to resemble this more:

but I'm not sure what is wrong with the way that I added the contours so that it's not doing that..

[SorooshMani-NOAA]

I think if you increase the expansion rates, you should be able to get closer to the last figure in terms of sizing. For the size function, you'll still get nan-data in the nan-region, but that's OK, since that region is going to be excluded when you mesh the geometry using that size function. In general, to get stepwise expansion moving away from a shoreline, I'd use constant value along with contour:

my_hfun = YOURHFUNTYPE(..., hmax=hmax, hmin=hmin)

# Just to use less memory for size function generation (so that you don't get small elements outside region of interest)
my_hfun.add_topo_bound_constraint(value=hmax, lower_bound=highest_level_for_meshing + eps, value_type='min')

my_hfun.add_constant_value(value=size_1, lower_bound=level_1)
my_hfun.add_contour(level=level_1, expansion_rate=0.05, target_size=size_1)
my_hfun.add_constant_value(value=size_2, lower_bound=level_2)
my_hfun.add_contour(level=level_2, expansion_rate=0.05, target_size=size_2)

The expansion_rate is something that you have to play with to get to the right value. If you want sharp transition, you can skip adding contours and just use constant value specifications.

For plotting the DEM, I don't have any specific recommendations, but I guess if you manage the nan/fill data in the contour's vmin/vmax you should get a better result.

[melissamattwig]

That makes sense, thank you!

I adjusted my code to include constant values along with the contours and changed the expansion rates. I figured it would take a while since it is a large DEM over Lake Ontario, but the code was stuck on adding the contour and constant value around the floodplain (the first contour and constant value added to the hfun after the topo bound constraints) for 18 hours when previously I was able to add 2 contours to the hfun in about 25 minutes. I looked at the raster.py source code to get an idea on how best to use the add_topo_bound_constraint argument but I'm not sure if that's conflicting with the contours or not. I thought that that argument sets limits for the element size based on the elevation level in the DEM (at least that's how I wrote it). Do you have an idea of what could be causing this to take so long when adding contours, constant values, and topo bound constraints or is 18+ hours to be expected when developing an hfun of this size? I'm running this code on our hpc.

print("adding contours")
print('adding topo bound constraint')
# Just to use less memory for size function generation (so that you don't get small elements outside region of interest)
hfun_r1.add_topo_bound_constraint(value=2700, lower_bound=ref_r1-245, value_type='max')
hfun_r1.add_topo_bound_constraint(value=35, upper_bound=ref_r1+60, value_type='min')
print('floodplain')
hfun_r1.add_constant_value(value=150, lower_bound=ref_r1+50)
hfun_r1.add_contour(level=ref_r1+50, expansion_rate=0.05, target_size=150)
print('shoreline')
hfun_r1.add_constant_value(value=100, lower_bound=ref_r1)
hfun_r1.add_contour(level=ref_r1, expansion_rate=0.05, target_size=100)
print('0m')
hfun_r1.add_constant_value(value=300, lower_bound=ref_r1-75)
hfun_r1.add_contour(level=ref_r1-75, expansion_rate=0.05, target_size=300)
print('75m depth')
hfun_r1.add_constant_value(value=800, lower_bound=ref_r1-150)
hfun_r1.add_contour(level=ref_r1-150, expansion_rate=0.05, target_size=800)
print('125m depth')
hfun_r1.add_constant_value(value=1500, lower_bound=ref_r1-200)
hfun_r1.add_contour(level=ref_r1-200, expansion_rate=0.05, target_size=1500)
print('150m depth')
hfun_r1.add_constant_value(value=2000, lower_bound=ref_r1-225)
hfun_r1.add_contour(level=ref_r1-225, expansion_rate=0.05, target_size=2000)
print('170m depth')
hfun_r1.add_constant_value(value=2700, lower_bound=ref_r1-245)
hfun_r1.add_contour(level=ref_r1-245, expansion_rate=0.05, target_size=2700)

[SorooshMani-NOAA]

hmm, that's a very long run time ... I know there are performance issues here and there but I'm not sure why it's taking that long. Looking at your script, I would first try to drop these:

hfun_r1.add_topo_bound_constraint(value=2700, lower_bound=ref_r1-245, value_type='max')
hfun_r1.add_topo_bound_constraint(value=35, upper_bound=ref_r1+60, value_type='min')

since they don't add much value as the hmax and hmin would do pretty much the same thing, however, you can instead add something like:

hfun_r1.add_topo_bound_constraint(value=2700, lower_bound=ref_r1+60, value_type='min')

or even +51 instead of +60 if your domain of interest ends at 50, so that you ensure all the nan value regions get a large mesh in the hfun (since hfun uses bbox).

The *constraint*s act as a barrier to the element size (instead of enforcing a size), so what you do with those two lines is to say:

1. any element above the deepest level should not be larger than 2700 m
2. any element below the shallowest level should not be smaller than 35 m

which is the same as setting hmin and hmax, but what we would like to achieve (for memory efficiency) is to say:

• Anything above the shallowest level should not be smaller than 2700 m.

I think this should help resolve the memory issue and reduce the runtime significantly, if that's not the case I'd suggest that you share the DEM with me if possible for me to test on my side to see if I find any bugs.

Update
Actually let's try with 2 contour+constants and then increase to the full script, just to see the timing instead of maybe waiting another 18+ hours before finding out! So I'd suggest replace the constraint lines as I suggested, then just keep two of the contour+constant levels and run and see how it goes, then if it goes through quickly try the full script.

[melissamattwig]

Unfortunately I ran the code again with just 2 contours/constants and so far it's been stuck in the same spot as before for about 2.5 hours. I tried it again with 2 different contours/constants to see if that will work but so far it seems like it may also take a long time.

What is the easiest way to share the DEM with you? It's about 2.7 gb

[SorooshMani-NOAA]

You can share it via Google drive or something like that if you have access, you can then email the link to me

[SorooshMani-NOAA]

@melissamattwig I tried running the script you provided and ran into memory issue (on my laptop). As you said the DEM is very high resolution (~10k x 40k pixels), which also results in very slow operations in general. The first thing I tried was to try processing the raster in chunks (passing chunk_size=5000 and overlap=100 to Raster constructor):

r1 = ocsmesh.Raster(
    'from_Melissa/LakeOntarioSeamlessDEMMeshTeamwithStLaurencev3PE15m.tif',
    chunk_size=5000,
    overlap=100)

It successfully went through

reading in USACE grid for domain
reading in raster
creating size function from raster
adding contours
adding topo bound constraint
floodplain
shoreline                          
adding shoreline constant

but failed at

adding shoreline contour

with a weird error related to multiprocessing (I'm still not sure why). I suggest that you first try the chunking as described above - maybe even try chunks of size 2500 - and see if you have success. If not, a quick and dirty way would be to reduce the DEM resolution to generate the size function, and then after you have the mesh you can interpolate bathy values from the original (hires) DEM.

To reduce the DEM resolution in OCSMesh you can use method:

OCSMesh/ocsmesh/raster.py

Lines 1153 to 1156 in 0563a79

	def resample(self,
	scaling_factor: float,
	resampling_method: Optional[str] = None
	) -> None:

I'll wait until you try these out, if these don't help we can again think of other solutions. Right now I don't have a better solution in mind.

[SorooshMani-NOAA]

Get bugfix/minor branch please

[melissamattwig]

I tried running the code with the chunking after updating ocsmesh with the chunking fix and I am also getting an error when trying to put in contours, except my error seems to be when I try to add the first contour and your error is with the second contours (after adding the floodplain contour). The error I'm getting is related to an empty linestring
raise EmptyPartError( shapely.errors.EmptyPartError: Can't create MultiLineString with empty component which I never got before the chunking. Do you think this is a bug related to the chunking?

[SorooshMani-NOAA]

I don't like the fact that we're getting different errors!! I might need to run the whole meshing process on my side once, when I have available resources. I'd suggest that you try the resample approach just to get some results. You can use a resampling factor such that the DEM resolution is still higher than your hmin or hmin/2. Let's see how that goes.

I created a ticket for myself (#130) to test the chunking approach on your DEM, but given my current available resources, I cannot promise any dates for running it.

[melissamattwig]

Thank you!! I also have a tangential question that I can pose in another thread if you'd prefer, but I tried to mesh lidar DEM files around Lake Ontario and received an error message about having too many file:
File "rasterio/_io.pyx", line 1491, in rasterio._io.DatasetWriterBase.__init__ rasterio.errors.RasterioIOError: Attempt to create new tiff file '/tmp/ocsmesh/5hnnf5kn' failed: Too many open files

I was trying to read in about 1300 1m lidar DEM files to the size function. I'm wondering if you've run into this issue before?

[SorooshMani-NOAA]

This should be related to ulimit set on your environment (if it's Linux). If you type ulimit -n, it should tell you how many files can be simultaneously opened. To modify that, just pass a new number, e.g. ulimit -n 2000. Then run the Python code in the same terminal session where you set the new ulimit for number of files.

[SorooshMani-NOAA]

@melissamattwig I was able to get a size function like the snapshot in less than an hour on my laptop.

I added a feature in OCSMesh that you need to add in your meshing script to generate the size function faster. I'm still testing it, when it's ready I'll let you know here.

[SorooshMani-NOAA]

I merged some bugfixes + performance improvements to the main branch today. You should be able to get an hfun like above in ~40min or so (at least on my laptop that's the timing). To get this performance, there are a couple of changes you need to make in your script. First of all for adding all the constraints and mesh size controls, you need to use the newly added hold_applying_added_constratins "Context Manger":

with hfun_r1.hold_applying_added_constratins():
    print("adding contours")
    print('adding topo bound constraint')
    hfun_r1.add_topo_bound_constraint(
        value=2700, lower_bound=ref_r1+51, value_type='min')
    print('floodplain')
    hfun_r1.add_constant_value(value=150, lower_bound=ref_r1+50)
    hfun_r1.add_contour(level=ref_r1+50, expansion_rate=0.05, target_size=150)
    print('shoreline')
    print('adding shoreline constant')
    hfun_r1.add_constant_value(value=100, lower_bound=ref_r1)
    # (profile1) First profiling stop
    print('adding shoreline contour')
    hfun_r1.add_contour(level=ref_r1, expansion_rate=0.05, target_size=100)
    print('0m')
    hfun_r1.add_constant_value(value=300, lower_bound=ref_r1-75)
    hfun_r1.add_contour(level=ref_r1-75, expansion_rate=0.05, target_size=300)
    print('75m depth')
    hfun_r1.add_constant_value(value=800, lower_bound=ref_r1-150)
    hfun_r1.add_contour(level=ref_r1-150, expansion_rate=0.05, target_size=800)
    print('125m depth')
    hfun_r1.add_constant_value(value=1500, lower_bound=ref_r1-200)
    hfun_r1.add_contour(level=ref_r1-200, expansion_rate=0.05, target_size=1500)
    print('150m depth')
    hfun_r1.add_constant_value(value=2000, lower_bound=ref_r1-225)
    hfun_r1.add_contour(level=ref_r1-225, expansion_rate=0.05, target_size=2000)
    print('170m depth')
    hfun_r1.add_constant_value(value=2700, lower_bound=ref_r1-245)
    hfun_r1.add_contour(level=ref_r1-245, expansion_rate=0.05, target_size=2700)

Note that the use of this hold method is crucial for performance (so that it doesn't reapply the constraints over and over again!

Another point is when you want to add any controls to the size function, you need to make sure matplotlib library which is used for extracting the contours from rasters is using agg (non interactive) backend. To do so, at the top of your script add:

import matplotlib
matplotlib.use('agg')

This is not needed all the times, but if you don't add it, then sometimes it might fail to get the contour refinement properly. Also note that when you add this you might not be able to create plots like before, so I'd suggest separating any plotting functions from the meshing or size function creation parts of the code. After done with mesh gen, you can then reset the backend to an interactive one.

In my case I just executed the script and then saved the mesh to a 2dm file and viewed on QGIS.
Please let me know how it goes.

See #132 #131 #130

[melissamattwig]

Thank you for your help with this!! I'm still getting the "can't create multiline string with empty component and I'm not sure why. I updated ocsmesh main branch using pip install and adding both code snippets for the hold and the matplotlib package. Here is what the error code says:

Traceback (most recent call last):
  File "/mnt/projects/hpc/mattwig/grid/ocsmesh/usace_dem.py", line 123, in <module>
    hfun_r1.add_contour(level=ref_r1+50, expansion_rate=0.05, target_size=150)
  File "/home/mattwig/miniconda3/envs/ocsmesh/lib/python3.11/site-packages/ocsmesh/features/constraint.py", line 403, in wrapped
    rv = method(obj, *args, **kwargs)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mattwig/miniconda3/envs/ocsmesh/lib/python3.11/site-packages/ocsmesh/hfun/raster.py", line 1022, in add_contour
    _contours = self.raster.get_contour(_level)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mattwig/miniconda3/envs/ocsmesh/lib/python3.11/site-packages/ocsmesh/raster.py", line 1371, in get_contour
    return self._get_raster_contour_feathered(level, iter_windows)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mattwig/miniconda3/envs/ocsmesh/lib/python3.11/site-packages/ocsmesh/raster.py", line 1497, in _get_raster_contour_feathered
    results = self._get_raster_contour_feathered_internal(
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mattwig/miniconda3/envs/ocsmesh/lib/python3.11/site-packages/ocsmesh/raster.py", line 1559, in _get_raster_contour_feathered_internal
    features = ops.linemerge(features)
               ^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mattwig/miniconda3/envs/ocsmesh/lib/python3.11/site-packages/shapely/ops.py", line 108, in linemerge
    source = MultiLineString([ls.coords for ls in lines])
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mattwig/miniconda3/envs/ocsmesh/lib/python3.11/site-packages/shapely/geometry/multilinestring.py", line 52, in __new__
    raise EmptyPartError(
shapely.errors.EmptyPartError: Can't create MultiLineString with empty component

[SorooshMani-NOAA]

That's strange ... what chunking parameters are you using? I mean for chunk_size and overlap? I'm using

    r1 = ocsmesh.Raster(
        'LakeOntarioSeamlessDEMMeshTeamwithStLaurencev3PE15m.tif',
        chunk_size=2500,
        overlap=100)

If you're using different values, I suggest you first use my values and see if it works, otherwise I can potentially add some workaround in a bugfix branch so that it skips empty parts to see if the rest of the process works for you or not.

[melissamattwig]

I'm using the same values as you... Here is my code up through adding contours and such:

## read in existing USACE mesh for domain
print("reading in USACE grid for domain")
mesh_domain = ocsmesh.Mesh.open('o4-v01-0m_igld85_ocsmesh.2dm', crs=4326)
## convert mesh domain into geom object
geom_domain = ocsmesh.Geom(mesh_domain)
## read in USACE lake ontario DEM
print("reading in raster")
#r1 = ocsmesh.Raster('data/LakeOntarioSeamlessDEMMeshTeamwithStLaurencev3PE15m.tif')
r1 = ocsmesh.Raster(
    'data/LakeOntarioSeamlessDEMMeshTeamwithStLaurencev3PE15m.tif',
    chunk_size=2500,
    overlap=100)
# size functon for just dem
print("creating size function from raster")
hfun_r1 = ocsmesh.Hfun(r1,hmin=30,hmax=3500)
# size function for multiple dems by creating list of rasters and passing through hfun
#hfun_r1 = ocsmesh.Hfun([other_rasters,r1],hmin=40,hmax=2500)

# set value for mean surface level for Lake Ontario
ref_r1 = 74.2 ##coastline

with hfun_r1.hold_applying_added_constratins():
    print("adding contours")
    print('adding topo bound constraint')
    hfun_r1.add_topo_bound_constraint(
        value=2700, lower_bound=ref_r1+51, value_type='min')
    print('floodplain')
    hfun_r1.add_constant_value(value=150, lower_bound=ref_r1+50)
    hfun_r1.add_contour(level=ref_r1+50, expansion_rate=0.05, target_size=150)
    print('shoreline')
    print('adding shoreline constant')
    hfun_r1.add_constant_value(value=100, lower_bound=ref_r1)
    # (profile1) First profiling stop
    print('adding shoreline contour')
    hfun_r1.add_contour(level=ref_r1, expansion_rate=0.05, target_size=100)
    print('0m')
    hfun_r1.add_constant_value(value=300, lower_bound=ref_r1-75)
    hfun_r1.add_contour(level=ref_r1-75, expansion_rate=0.05, target_size=300)
    print('75m depth')
    hfun_r1.add_constant_value(value=800, lower_bound=ref_r1-150)
    hfun_r1.add_contour(level=ref_r1-150, expansion_rate=0.05, target_size=800)
    print('125m depth')
    hfun_r1.add_constant_value(value=1500, lower_bound=ref_r1-200)
    hfun_r1.add_contour(level=ref_r1-200, expansion_rate=0.05, target_size=1500)
    print('150m depth')
    hfun_r1.add_constant_value(value=2000, lower_bound=ref_r1-225)
    hfun_r1.add_contour(level=ref_r1-225, expansion_rate=0.05, target_size=2000)
    print('170m depth')
    hfun_r1.add_constant_value(value=2700, lower_bound=ref_r1-245)
    hfun_r1.add_contour(level=ref_r1-245, expansion_rate=0.05, target_size=2700)

[SorooshMani-NOAA]

OK, I'll let you know once I've prepared the workaround on a bugfix branch