Error when running mesh.interpolate(raster)

[TPCollings]

Hi there,

I get an error when running the mesh.interpolate(raster) function. I'm following an example in the documentation to create a base mesh before further refinement. My code is below and the error are included below.

Thanks for your help,

Tom

import pathlib
import os 
from copy import deepcopy
import numpy as np
import geopandas as gpd
from shapely import geometry
import jigsawpy

from ocsmesh import Raster, Geom, Hfun, Mesh, JigsawDriver, utils

dempath = pathlib.Path('/home/tom/fathom-tide-surge-model/data/bathy/GEBCO_2022.nc')
domainpath = '/home/tom/fathom-tide-surge-model/data/domains/N_Atl01_shp/N_Atl_boundary_01.shp'

domain = gpd.read_file(domainpath)

raster1 = Raster(dempath) 
raster2 = Raster(dempath)
raster3 = Raster(dempath)

base_geom = Geom([raster1],
                 base_shape=domain.unary_union,
                 base_shape_crs=domain.crs,
                 zmax=10)

base_hfun = Hfun([raster2],
                 base_shape=domain.unary_union,
                 base_shape_crs=domain.crs,
                 hmin=1000,hmax=15000,
                 method='fast')

driver = JigsawDriver(geom=base_geom,hfun=base_hfun)
base_mesh=driver.run()
base_mesh.interpolate(raster3)

---------------------------------------------------------------------------
MaybeEncodingError                        Traceback (most recent call last)
Input In [97], in <cell line: 1>()
----> 1 base_mesh.interpolate(raster3)

File ~/anaconda3/envs/ocsmesh_main1/lib/python3.9/site-packages/ocsmesh/mesh/mesh.py:369, in EuclideanMesh2D.interpolate(self, raster, method, nprocs)
    367 if nprocs > 1:
    368     with Pool(processes=nprocs) as pool:
--> 369         res = pool.starmap(
    370             _mesh_interpolate_worker,
    371             [(self.vert2['coord'], self.crs,
    372                 _raster.tmpfile, _raster.chunk_size, method)
    373              for _raster in raster]
    374             )
    375     pool.join()
    376 else:

File ~/anaconda3/envs/ocsmesh_main1/lib/python3.9/multiprocessing/pool.py:372, in Pool.starmap(self, func, iterable, chunksize)
    366 def starmap(self, func, iterable, chunksize=None):
    367     '''
    368     Like `map()` method but the elements of the `iterable` are expected to
    369     be iterables as well and will be unpacked as arguments. Hence
    370     `func` and (a, b) becomes func(a, b).
    371     '''
--> 372     return self._map_async(func, iterable, starmapstar, chunksize).get()

File ~/anaconda3/envs/ocsmesh_main1/lib/python3.9/multiprocessing/pool.py:771, in ApplyResult.get(self, timeout)
    769     return self._value
    770 else:
--> 771     raise self._value

MaybeEncodingError: Error sending result: '<multiprocessing.pool.ExceptionWithTraceback object at 0x7fca6c612460>'. Reason: 'PicklingError("Can't pickle <class 'dfitpack.error'>: import of module 'dfitpack' failed")'

[SorooshMani-NOAA]

Hi @TPCollings, let me run a simple interpolation test case and check if I see the same issue. If not, I need the shapefile you're using to run exactly what you're running to see why you get this error.

Without any testing, it seems to me that something unexpected happens in the internal interpolation code (at scipy level). This still could be from how OCSMesh is handling the arguments passed to interpolation code or something that is causing problem at the topobathy file. Would it be possible that you share the shapefile you're using with me?

[TPCollings]

Hi @SorooshMani-NOAA, thanks for your reply. Yeah I've attached the shapefile here.

Thanks,

Tom N_Atl01_shp.zip

[SorooshMani-NOAA]

@TPCollings thanks for sharing the file. May I ask if you installed OCSMesh from Github repo or from PyPI? We haven't release the latest Github repo on PyPI yet, so there might be fixes that didn't make it to the PyPI version

[TPCollings]

@SorooshMani-NOAA Ah, I installed via PyPi. I'll try again using the Github repo and see if that helps.

[TPCollings]

@SorooshMani-NOAA Before downloading the repo and installing with the commands given in the ReadMe on Github, should I set up the conda environment using the .yml file as before?

[SorooshMani-NOAA]

Your test case is still running on my machine, so I don't know if getting latest is necessarily going to fix the issue. But to answer your question, yes, you need to have your conda environment created before installing OCSMesh from the repo as well. Getting all the libraries to play nice without conda can be a headache. I would also suggest that you use mamba package to resolve the conda environment faster. You can install mamba in your conda base environment by:

conda install -n base -cconda-forge mamba

and then later you can just use mamba to run install or create environment commands (I still don't use it for activate or deactivate commands of the environment). So you would do something like:

mamba env create -f environment.yml

[SorooshMani-NOAA]

@TPCollings I got the same error you have. Let me dig a little deeper and see what is going on

[TPCollings]

Okay, thank you very much!

[SorooshMani-NOAA]

Running the test case with base_mesh.interpolate(raster3, nprocs=1) I get an error related to the interpolation function in scipy:

Traceback (most recent call last):
  File "/home/Soroosh.Mani/test_ocsmesh.py", line 36, in <module>
    base_mesh.interpolate(raster3, nprocs=1)
  File "/home/Soroosh.Mani/sandbox/ocsmesh/ocsmesh/mesh/mesh.py", line 385, in interpolate
    res = [_mesh_interpolate_worker(
  File "/home/Soroosh.Mani/sandbox/ocsmesh/ocsmesh/mesh/mesh.py", line 385, in <listcomp>
    res = [_mesh_interpolate_worker(
  File "/home/Soroosh.Mani/sandbox/ocsmesh/ocsmesh/mesh/mesh.py", line 2079, in _mesh_interpolate_worker
    f = RectBivariateSpline(
  File "/home/Soroosh.Mani/miniconda3/envs/ocsmesh/lib/python3.9/site-packages/scipy/interpolate/fitpack2.py", line 1345, in __init__
    nx, tx, ny, ty, c, fp, ier = dfitpack.regrid_smth(x, y, z, xb, xe, yb,
dfitpack.error: (len(z)==mx*my) failed for 3rd argument z

I need to check and see why this happens ...

[SorooshMani-NOAA]

Might be related to https://github.com/scipy/scipy/issues/7663

[SorooshMani-NOAA]

@TPCollings, can you please try this on a smaller domain size and see if the issue is resolved or not? If that's the issue then you can specify a chunk_size when you create a raster object for raster3 so that raster is interpolated piece by piece.

Please let me know if this works for your case, thanks!

[TPCollings]

@SorooshMani-NOAA simply reducing the domain size hasn't worked. I tried using a 2 different sizes, one that covered the Gulf of Mexico, and one that covered New Orleans, but neither worked.

I'm trying to use the chunk_size flag to see if it works with that. What are the units of the chunk_size? What would be an appropriate value? Thanks

[SorooshMani-NOAA]

@TPCollings the chunk_size is just the number of raster pixels to consider for windowing/chunking of the DEM; it results in creating square window chunks from the raster to read it piece by piece.

The fact that smaller size doesn't work however, makes me think it might be something else ... can you try interpolating with .interpolate(raster3, nprocs=1) and see what error you get? Is it the same thing I see in this comment?

[SorooshMani-NOAA]

@TPCollings I forgot to ask, did you try other interpolation methods as well? e.g. linear or nearest? You can pass the interpolation method to the .interpolate(...) call as an argument e.g. method='linear'. spline is the default.

[TPCollings]

@SorooshMani-NOAA So I've switched to using a .tif version of the GEBCO dem, rather than .nc, and this has worked. I haven't tried it on a large domain, but it worked for a small test area around New Orleans.

[SorooshMani-NOAA]

Thanks for the update. That's curious. Maybe the netcdf reader of rasterio causes some issue somewhere. The code in ocsmesh doesn't really care whether it's nc or tif. Although most of the test cases run using GeoTiff.