cnes / pangeo-pyinterp

Python library for optimized interpolation.

pangeo-pyinterp

Python library for optimized geo-referenced interpolation.

About

The motivation of this project is to provide tools for interpolating geo-referenced data used in the field of geosciences. Other libraries cover this problem, but written entirely in Python, the performance of these projects was not quite sufficient for our needs. That is why this project started.

With this library, you can interpolate 2D, 3D, or 4D fields using n-variate and bicubic interpolators and unstructured grids. You can also apply for a data binning on the bivariate area by simple or linear binning.

The library core is written in C++ using the Boost C++ Libraries, Eigen3 and pybind11 libraries.

This software also uses CMake to configure the project and Googletest to perform unit testing of the library kernel.

Fill undefined values

The undefined values in the grids do not allow interpolation of values located in the neighborhood. This behavior is a concern when you need to interpolate values near the mask of some fields. The library provides utilities to fill the undefined values:

• loess to fill the undefined values on the boundary between the defined/undefined values using local regression.
• gauss_seidel to fill all undefined values in a grid using the Gauss-Seidel method by relaxation.

Geographic indexers

N-Dimensional Grids

N-dimensional grid is a grid defined by a matrix, in a 2D space, by a cube in a 3D space, etc. Each dimension of the grid is associated with a vector corresponding to its coordinates or axes. Axes used to locate a pixel in the grid from the coordinates of a point. These axes are either:

• regular: a vector of 181 latitudes spaced a degree from -90 to 90 degrees;
• irregular: a vector of 109 latitudes irregularly spaced from -90 to 89.940374 degrees.

These objects are manipulated by the class pyinterp.Axis, which will choose, according to Axis definition, the best implementation. This object will allow you to find the two indexes framing a given value. This operating mode allows better performance when searching for a regular axis (a simple calculation will enable you to see the index of a point immediately). In contrast, in the case of an irregular axis, the search will be performed using a binary search.

Finally, this class can define a circular axis from a vector to correctly locate a value on the circle. This type of Axis will is used for handling longitudes.

Temporal Axes

The pyinterp.TemporalAxis class handles temporal axes, i.e., axes defined by 64-bit integer vectors, which is the encoding used by numpy to control dates. This class allows you to process dates using integer arithmetic to ensure that no information is lost during calculations. These objects are used by spatiotemporal grids to perform temporal interpolations.

Unstructured Grids

In the case of unstructured grids, the index used is a R*Tree. These trees have better performance than the KDTree generally found in Python library implementations.

The tree used here is the implementation provided by the C++ Boost library.

An adaptation has been introduced to address spherical equatorial coordinates effectively. Although the Boost library allows these coordinates to manipulate natively, the performance is lower than in the case of Cartesian space. Thus, we have chosen to implement a conversion of Longitude Latitude Altitude (LLA) coordinates into Earth-Centered, Earth-Fixed (ECEF) coordinates transparently for the user to ensure that we can preserve excellent performance. The disadvantage of this implementation is that it requires fairly more memory, as one more element gets used to index the value of the Cartesian space.

The management of the LLA/ECEF coordinate conversion is managed to use the Olson, D.K. algorithm. It has an excellent performance with the accuracy of 1e-8 meters for altitude.

Geohash

Geohashing is a geocoding method used to encode geographic coordinates (latitude and longitude) into a short string of digits and letters delineating an area on a map, which is called a cell, with varying resolutions. The more characters in the string, the more precise the location.

Geohashes use Base-32 alphabet encoding (characters can be 0 to 9 and A to Z, excl A, I, L and O).

The geohash is a compact way of representing a location, and is useful for storing a location in a database, or for indexing a location in a database.

pangeo-pyinterp's Issues

Extract lat, lon from RTree

Is your feature request related to a problem? Please describe.
When using RTree, we can extract distances and values of nearest neighbors (for example when using the query method). Currently, there is no way to get the (lon, lat) position of the nearest neighbors through the method.

I am interested in finding the (lon, lat) coordinates of the nearest neighbors using any of the RTree methods.

Describe the solution you'd like
Return the position in the query method in addition to the standard output of distance in meters and the value of the different neighbors.

Describe alternatives you've considered
Loop through the data without RTree to calculate the spheroidal distances to the point of interest ($$$), then compare distances with query method results to save (lon, lat) pairs 💀

Additional context
There may already be a process for this, but I couldn't find anything in the PyInterp Read the Docs

Creating a 3D interpolator from a large dataset results in a much larger ram spike

Describe the bug
I want to create an interpolator from a large number of netcdf files. The dataset size is 10.3GB, however it contains 10 variables so each one takes 5.1GB. (I have 62GB of ram available so pretty comfortably fits within memory). When I go to create the Grid3D object, the memory usage spikes and I quickly run out of memory and the program stops.

Is there advise on how much more memory is used when creating the Grid3D object, compared to the datset size?

Do you know of any optimisations/tricks in xarray to help this? (I tried playing with chunk sizes but I couldnt get anything on that end to work)

Does the ram usage drop back down to roughly 10GB after the object has been created?

To Reproduce

import pyinterp
import xarray

ds = xarray.open_mfdataset(files, parallel=True)
ds = ds.load()  # not required, but just proving that it can be fit in memory
ds["lon"] = ds.lon.assign_attrs(
    {"standard_name": "longitude", "long_name": "longitude", "units": "degrees_east", "axis": "X"}
)
ds["lat"] = ds.lat.assign_attrs(
    {"standard_name": "latitude", "long_name": "latitude", "units": "degrees_north", "axis": "Y"}
)
ds["time"] = ds.time.assign_attrs({"standard_name": "time", "axis": "T"})

If required, I can upload the netcdf files to fully replicate

Expected behavior
Expected behavior is memory usage to increase, but hopefully still be able to create the object with 64 gigs of ram

Pyinterp/Numpy/Python version information
Output from

Pyinterp version: 2022.9.1
NumPy version: 1.24.2
Python version: 3.10.8

Support for Apple silicon

Hope to support apple m1.

Feature request: fill.loess that accepts Grid4D

Is your feature request related to a problem? Please describe.

Hello, first of, great tool. I have been using pyinterp for some time to interpolate sea surface height data (dims: time/lat/lon) without any problems. But now that I need to use it with ocean current data (dims: time/depth/lat/lon) I am having some difficulties. To avoid creating NaNs values near land, I always fill some NaNs before the interpolation. That works fine with 2D and 3D grids, but with a 4D grid I get the following error:

da = ds['sw_cur_u']
grid = pyinterp.backends.xarray.Grid4D(da)
filled = pyinterp.fill.loess(grid, nx=3, ny=3)
...
TypeError: loess_float32(): incompatible function arguments. The following argument types are supported:
    1. (grid: pyinterp.core.Grid2DFloat32, nx: int = 3, ny: int = 3, value_type: pyinterp.core.fill.ValueType = <ValueType.Undefined: 0>, num_threads: int = 0) -> numpy.ndarray[numpy.float32]
    2. (grid: pyinterp.core.Grid3DFloat32, nx: int = 3, ny: int = 3, value_type: pyinterp.core.fill.ValueType = <ValueType.Undefined: 0>, num_threads: int = 0) -> numpy.ndarray[numpy.float32]
    3. (grid: pyinterp.core.TemporalGrid3DFloat32, nx: int = 3, ny: int = 3, value_type: pyinterp.core.fill.ValueType = <ValueType.Undefined: 0>, num_threads: int = 0) -> numpy.ndarray[numpy.float32]

Invoked with: <pyinterp.core.TemporalGrid4DFloat32 object at 0x7ff5a6c4cb30>, 3, 3, <ValueType.Undefined: 0>, 0

Describe the solution you'd like

A pyinterp.fill.loess that accepts Grid4D data. Only filling along the x (longitude) and y (latitude) dims is enough, no need to also fill along z (depth/height) and t (time) dimensions.

Describe alternatives you've considered

Iterating over the depth dimension and joining all the data arrays after. It is possible, but a little convoluted, specially when trying to write a flexible function that accepts 2D, 3D and 4D xarray dataarrays (I am trying to do just that), e.g.:

for var in ds_old.data_vars:
    grid = myclass(grid_old, grid_new)
    ds_new[var] = grid.interp(ds_old[var], fill_nans=4)

Additional context

Here an example of a netcdf file with 4D variables (time,depth,lat,lon) that fails when filling NaNs.
hycom.zip

Thank you.

Python kernel crashes when running "ex_geohash.py" example

Describe the bug
Python kernel crashes when running ex_geohash example..

python(23293,0x1159de600) malloc: *** error for object 0x7ff7b6c160a8: pointer being freed was not allocated
python(23293,0x1159de600) malloc: *** set a breakpoint in malloc_error_break to debug
zsh: abort      python ex_geohash.py

Pyinterp/Numpy/Python version information
I run Python 3.10.6 with the following packages

  - cartopy==0.21.0
  - contourpy==1.0.5
  - cycler==0.11.0
  - fonttools==4.38.0
  - gdal==3.4.1 and libgdal=3.4.1
  - kiwisolver==1.4.4
  - matplotlib==3.6.1
  - pillow==9.2.0
  - pyproj==3.4.0
  - pyqt5-sip==12.11.0
  - pyshp==2.3.1
  - shapely==1.8.5.post1
  - numpy=1.23.4
  - ipykernel=6.15.2
  - ipython=8.4.0

Additional context
Add any other context about the problem here.

within functionality flipped for IDW

Describe the bug
within kwarg seems to be flipped for inverse_distance_weighting and in example notebooks. See

pangeo-pyinterp/src/pyinterp/core/include/pyinterp/detail/geometry/rtree.hpp

Line 184 in 47520d6

auto nearest = within ? query(point, k) : query_within(point, k);

Also, in https://github.com/CNES/pangeo-pyinterp/blob/master/notebooks/auto_examples/ex_unstructured.ipynb there are comments saying that within = False implies "# Extrapolation is forbidden". In reality, this should imply that extrapolation is allowed.

Same with https://github.com/CNES/pangeo-pyinterp/blob/master/notebooks/auto_examples/pangeo_unstructured_grid.ipynb

Expected behavior
within=True should imply that extrapolation is forbidden, which matches the docstrings for the functions and matches the behavior for radial_basis_function and query

Pyinterp/Numpy/Python version information
Output from

import sys, numpy, pyinterp.version
print(pyinterp.version.release(True))
print(numpy.__version__)
print(sys.version)

I get an AttributeError b/c pyinterp has no attribute version. But pyinterp.__version__ yields 0.6.1

Provide manylinux wheels

Is your feature request related to a problem? Please describe.
Providing pre-built wheels built inside a manylinux container makes the package easier and faster to install on Linux.

Describe the solution you'd like
Create python wheels inside a manylinux container and publish those to pypi.

Describe alternatives you've considered
I could build the package myself inside a manylinux container and move it onto my server but may as well add it to the main package if possible :)

Additional context
I have done this sort of thing before but I haven't ran it with github workflows (i've done it with Jenkins).

There are some example workflows here:

• https://github.com/scikit-image/scikit-image/blob/main/.github/workflows/wheel_tests_and_release.yml
• https://github.com/scipy/scipy/blob/main/.github/workflows/wheels.yml

And an example in the cibuildwheel package

name: Build and upload to PyPI

# Build on every branch push, tag push, and pull request change:
on: [push, pull_request]
# Alternatively, to publish when a (published) GitHub Release is created, use the following:
# on:
#   push:
#   pull_request:
#   release:
#     types:
#       - published

jobs:
  build_wheels:
    name: Build wheels on ${{ matrix.os }}
    runs-on: ${{ matrix.os }}
    strategy:
      matrix:
        os: [ubuntu-20.04, windows-2019, macos-11]

    steps:
      - uses: actions/checkout@v3

      - name: Build wheels
        uses: pypa/[email protected]

      - uses: actions/upload-artifact@v3
        with:
          path: ./wheelhouse/*.whl

  build_sdist:
    name: Build source distribution
    runs-on: ubuntu-latest
    steps:
      - uses: actions/checkout@v3

      - name: Build sdist
        run: pipx run build --sdist

      - uses: actions/upload-artifact@v3
        with:
          path: dist/*.tar.gz

  upload_pypi:
    needs: [build_wheels, build_sdist]
    runs-on: ubuntu-latest
    # upload to PyPI on every tag starting with 'v'
    if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/v')
    # alternatively, to publish when a GitHub Release is created, use the following rule:
    # if: github.event_name == 'release' && github.event.action == 'published'
    steps:
      - uses: actions/download-artifact@v3
        with:
          # unpacks default artifact into dist/
          # if `name: artifact` is omitted, the action will create extra parent dir
          name: artifact
          path: dist

      - uses: pypa/[email protected]
        with:
          user: __token__
          password: ${{ secrets.pypi_password }}
          # To test: repository_url: https://test.pypi.org/legacy/

I am happy to try implement this myself and do a PR.

Understanding trivariate interpolator

I‘m looking for something similar to the regular grid interpolator in SciPy [1]. How does the trivariate interpolator compare, algorithmically and from practical application perspective?

[1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html#scipy.interpolate.RegularGridInterpolator

Can't import pyinter.fill

I'm working on Windows 10 and I've installed the Boost C++, Eigen3, GNU Scientific Library, and pybind11 libraries. Then I've proceeded to conda install on the Anacona Prompt the pyinterp library. But when I intended to import some functions, the next occurred:

import pyinterp.fill

This is the complete traceback:
<<
Cell In[2], line 1
import pyinterp.fill

File ~\anaconda3\lib\site-packages\pyinterp_init_.py:9
from . import geodetic, geohash, version

File ~\anaconda3\lib\site-packages\pyinterp\geodetic_init_.py:9
from .. import interface

File ~\anaconda3\lib\site-packages\pyinterp\interface.py:14
from . import core

ImportError: DLL load failed while importing core: the specified module can not be found.

I've searched online about what DLL file could be missing, but no responses yet.

I would appreciate if you could guide me about what files I might be missing or if it is necessary to uninstall it and reinstall it in a different way, thank you very much in advance

Using Euclidean Coordinates

Hi- your interpolation library is one of the fastest around with the C++ implementation. I was hoping to use some of these functions for a project where I have a large shape=(1540, 1049, 521) array with scattered data of 3D vectors (~10,000 scattered points). Your functions work very quickly, but the casting to spherical coordinates is throwing off my results! Is there a manual way around this? I see the system parameter for Rtree, but I'm having a hard time figuring out what the valid options are.
Thanks for your time!

Here is my code snippet- force_locations is an mx3 array, where m is the number of scattered data

import pyinterp
field = np.zeros((vessel_seg.shape[0], vessel_seg.shape[1], vessel_seg.shape[2], 3))
n, p, q = vessel_seg.shape

print("Interpolating IDW")
mesh_x = pyinterp.RTree()

mesh_x.packing(force_locations[:], force_displacements[:, 0])

mx, my, mz = np.meshgrid(np.arange(n),
np.arange(p),
np.arange(q))
# indexing='ij')

coords = np.vstack((mx.ravel(), my.ravel(), mz.ravel())).T

rbf_x, neighbors = mesh_x.radial_basis_function(
coords,
within=False,
k=19,
rbf='cubic',
# smooth=1e-4,
num_threads=0)
field[coords[:, 0], coords[:, 1], coords[:, 2], 0] = rbf_x

Error in pyinterp.core when trying to import pyinterp on Windows

Describe the bug
When trying to import pyinterp (v2023.2.1), I get the following error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\manow\miniconda3\envs\py38_ml_TfGpuMain\lib\site-packages\pyinterp\__init__.py", line 9, in <module>
    from . import geodetic, geohash, version
  File "C:\Users\manow\miniconda3\envs\py38_ml_TfGpuMain\lib\site-packages\pyinterp\geodetic\__init__.py", line 9, in <module>
    from .. import interface
  File "C:\Users\manow\miniconda3\envs\py38_ml_TfGpuMain\lib\site-packages\pyinterp\interface.py", line 14, in <module>
    from . import core
ImportError: DLL load failed while importing core: The specified module could not be found.

Looking at the core.cp38-win_amd64.pyd library with Dependencies, I see that the library is linked against gsl-25.dll. However, the gsl (v2.7.1) library is just gsl.dll

To Reproduce
Install pyinterp=2023.2.1=mkl_py38hb2a16ad_0 via conda on Windows 10
import pyinterp

Expected behavior
No crash

Pyinterp/Numpy/Python version information
Output from

import sys, numpy, pyinterp.version
print(pyinterp.__version__)
print(numpy.__version__)
print(sys.version)

pyinterp v2023.2.1
numpy v1.23.5
python v3.8.16 (default, Jan 17 2023, 22:25:28) [MSC v.1916 64 bit (AMD64)]

Wrong kwarg fitting_model value in bicubic function

Describe the bug

Hello, not a big problem, but I noticed that the docstring for the bicubic function shows bicubic as the default value for the fitting_model parameter (link), while in the function declaration, c_spline is the default (link).

I thought of doing a simple PR to fix this, but I am not sure which one is to be used as default.

Thank you.

Install CMakeLists.txt not found

Im running Ubuntu 18.04. I'm stuck on one error for install and hoping you can point me in the right direction. From the help files I can't seem to find a way to set a path, and/or can't figure our which library needs to be installed to add BLAS.

python setup.py build --boost-root /usr/lib/x86_64-linux-gnu/ --gsl-root /usr/lib/x86_64-linux-gnu/ --eigen-root /home/mmann1123/anaconda3/envs/geowombat/lib/python3.7/site-packages/scipy/sparse/linalg/eigen --mkl-root /home/mmann1123/anaconda3/envs/geowombat/lib/

running build
running build_py
copying src/pyinterp/version.py -> build/lib.linux-x86_64-3.7/pyinterp
running build_ext
cmake /home/mmann1123/Documents/github/pangeo-pyinterp -DCMAKE_BUILD_TYPE=Release -DCMAKE_LIBRARY_OUTPUT_DIRECTORY=/home/mmann1123/Documents/github/pangeo-pyinterp/build/lib.linux-x86_64-3.7/pyinterp -DPYTHON_EXECUTABLE=/home/mmann1123/anaconda3/envs/geowombat/bin/python -DBOOSTROOT=/usr/lib/x86_64-linux-gnu/ -DGSL_ROOT_DIR=/usr/lib/x86_64-linux-gnu/ -DEIGEN3_INCLUDE_DIR=/home/mmann1123/anaconda3/envs/geowombat/lib/python3.7/site-packages/scipy/sparse/linalg/eigen
-- Build type: Release
-- Boost version: 1.65.1
-- A library with BLAS API not found. Please specify library location.
CMake Error at CMakeLists.txt:190 (add_subdirectory):
The source directory

/home/mmann1123/Documents/github/pangeo-pyinterp/third_party/googletest

does not contain a CMakeLists.txt file.

CMake Error at CMakeLists.txt:200 (add_subdirectory):
The source directory

/home/mmann1123/Documents/github/pangeo-pyinterp/third_party/pybind11

does not contain a CMakeLists.txt file.

CMake Error at src/pyinterp/core/CMakeLists.txt:13 (pybind11_add_module):
Unknown CMake command "pybind11_add_module".

-- Configuring incomplete, errors occurred!
See also "/home/mmann1123/Documents/github/pangeo-pyinterp/build/temp.linux-x86_64-3.7/CMakeFiles/CMakeOutput.log".
See also "/home/mmann1123/Documents/github/pangeo-pyinterp/build/temp.linux-x86_64-3.7/CMakeFiles/CMakeError.log".
error: command 'cmake' failed with exit status 1

Interpolate for multiple data variables with a single interpolate call

Is your feature request related to a problem? Please describe.
I want to interpolate for multiple variables on the same grid (My grib data holds wind speed, wind direction U, wind direction V etc). Currently, I need to make an interpolator for each variable like:

interpolators = {
    'speed': pyinterp.backends.xarray.Grid3D(ds.speed),
    'u': pyinterp.backends.xarray.Grid3D(ds.u),
    'v': pyinterp.backends.xarray.Grid3D(ds.v),
}

So when I want to interpolate for given lat, lon, time I have to interpolate the grid three times even though they are at the same point in the 3D grid. I think it would be more efficient to interpolate each of the three variables simultaneously.

Describe the solution you'd like
Be able to make an interpolator from an XArray dataset, rather than a DataArray.

Example:

import numpy
import pyinterp.backends.xarray
import pyinterp.tests

ds = pyinterp.tests.load_grid3d()
ds['something_else'] = ds.time.astype(int)/1e16 + ds.latitude + ds.longitude
interpolator = pyinterp.backends.xarray.Grid3D(ds)

mx, my, mz = numpy.meshgrid(numpy.arange(-180, 180, 0.25) + 1 / 3.0,
                            numpy.arange(-80, 80, 0.25) + 1 / 3.0,
                            numpy.array(['2002-07-02T15:00:00'],
                                        dtype='datetime64'),
                            indexing='ij')

trivariate = interpolator.trivariate(
    dict(longitude=mx.ravel(), latitude=my.ravel(), time=mz.ravel()))

Which will result in something like:

{'tcw': array([ 1.00482318,  1.01251311,  1.02020304, ..., 11.12754145,
       11.10233118, 11.07064156]),
 'something_else': array([203.22886667, 203.47886667, 203.72886667, ..., 362.22886667,
       362.47886667, 362.72886667])
}

Describe alternatives you've considered
I tried to stack the data axis and use a 4d interpolator:

interpolator = pyinterp.backends.xarray.Grid4D(ds.to_array())

However, the variable axis is not supported (it is not np.float64)

I tried dropping the variable axis:

interpolator = pyinterp.backends.xarray.Grid4D(ds.to_array().drop("variable"))

Which makes the interpolator, but I cant seem to interpolate it anywhere:

interpolator.quadrivariate(dict(longitude=mx.ravel(), latitude=my.ravel(), time=mz.ravel()))

Traceback (most recent call last):
  File "/home/tom.mclean/miniconda3/envs/routingserver/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3378, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-38-859473664d87>", line 1, in <module>
    tmp.quadrivariate(dict(longitude=mx.ravel(), latitude=my.ravel(), time=mz.ravel()))
  File "/home/tom.mclean/miniconda3/envs/routingserver/lib/python3.10/site-packages/pyinterp/backends/xarray.py", line 390, in quadrivariate
    self, *_coords(coords, self._dims, self._datetime64), *args,
  File "/home/tom.mclean/miniconda3/envs/routingserver/lib/python3.10/site-packages/pyinterp/backends/xarray.py", line 129, in _coords
    raise IndexError('too many indices for array')
IndexError: too many indices for array

Or:

tmp.quadrivariate(dict(longitude=mx.ravel(), latitude=my.ravel(), time=mz.ravel(), variable="tcw"))
Traceback (most recent call last):
  File "/home/tom.mclean/miniconda3/envs/routingserver/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3378, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-39-3ee1e6de957d>", line 1, in <module>
    tmp.quadrivariate(dict(longitude=mx.ravel(), latitude=my.ravel(), time=mz.ravel(), variable="tcw"))
  File "/home/tom.mclean/miniconda3/envs/routingserver/lib/python3.10/site-packages/pyinterp/backends/xarray.py", line 389, in quadrivariate
    return interpolator.quadrivariate(
  File "/home/tom.mclean/miniconda3/envs/routingserver/lib/python3.10/site-packages/pyinterp/interpolator/quadrivariate.py", line 55, in quadrivariate
    return getattr(core, function)(instance,
TypeError: quadrivariate_float64(): incompatible function arguments. The following argument types are supported:
    1. (grid: pyinterp.core.Grid4DFloat64, x: numpy.ndarray[numpy.float64], y: numpy.ndarray[numpy.float64], z: numpy.ndarray[numpy.float64], u: numpy.ndarray[numpy.float64], interpolator: pyinterp.core.BivariateInterpolator3D, z_method: Optional[str] = None, u_method: Optional[str] = None, bounds_error: bool = False, num_threads: int = 0) -> numpy.ndarray[numpy.float64]
    2. (grid: pyinterp.core.TemporalGrid4DFloat64, x: numpy.ndarray[numpy.float64], y: numpy.ndarray[numpy.float64], z: numpy.ndarray[numpy.int64], u: numpy.ndarray[numpy.float64], interpolator: pyinterp.core.TemporalBivariateInterpolator3D, z_method: Optional[str] = None, u_method: Optional[str] = None, bounds_error: bool = False, num_threads: int = 0) -> numpy.ndarray[numpy.float64]
Invoked with: <pyinterp.core.TemporalGrid4DFloat64 object at 0x7feb47eb4d70>, array([-179.66666667, -179.66666667, -179.66666667, ...,  180.08333333,
        180.08333333,  180.08333333]), array([-79.66666667, -79.41666667, -79.16666667, ...,  79.58333333,
        79.83333333,  80.08333333]), array(['2002-07-02T15:00:00.000000000', '2002-07-02T15:00:00.000000000',
       '2002-07-02T15:00:00.000000000', ...,
       '2002-07-02T15:00:00.000000000', '2002-07-02T15:00:00.000000000',
       '2002-07-02T15:00:00.000000000'], dtype='datetime64[ns]'), array('tcw', dtype='<U3'), <pyinterp.core.TemporalBilinear3D object at 0x7feb8023b6f0>; kwargs: z_method='linear', u_method='linear', bounds_error=False, num_threads=0

Wrapping with a transpose matrix

Describe the bug
Wrong result with wrapping

To Reproduce

import pyinterp
import numpy as np

x_axis = np.arange(0,360,20)
y_axis = np.arange(-90,90,20)
z = np.outer(x_axis,y_axis).T
axis = [pyinterp.Axis(y_axis, is_circle=False), pyinterp.Axis(x_axis, is_circle=True)]
g = pyinterp.Grid2D(*axis, z)
x, y = np.array([20]), np.array([10])
print(z.shape)
print(x - 360, y, pyinterp.bivariate(g, y, x - 360, num_threads=1))
print(x, y, pyinterp.bivariate(g, y, x, num_threads=1))
print(x + 360, y, pyinterp.bivariate(g, y, x + 360, num_threads=1))

Output

(9, 18)
[-3400.]
[200.]
[3800.]

Expected behavior
Result must be 200 for each call.

Pyinterp/Numpy/Python version information

>>> print(pyinterp.__version__)
0.15.2
>>> print(numpy.__version__)
1.22.0
>>> print(sys.version)
3.10.4 (main, Mar 31 2022, 08:41:55) [GCC 7.5.0]
>>> 

ImportError: libomp.so: cannot open shared object file: No such file or directory

Hi:
I install pyinterp with "conda install pyinterp -c conda-forge" in ubuntu 19.10 system. When I import the library, I got the following errors:

>>> import pyinterp
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/xxx/anaconda3/lib/python3.7/site-packages/pyinterp/__init__.py", line 6, in <module>
    from .axis import TemporalAxis
  File "/home/xxx/anaconda3/lib/python3.7/site-packages/pyinterp/axis.py", line 12, in <module>
    from . import core
ImportError: libomp.so: cannot open shared object file: No such file or directory

So how can I solve this problem, thanks very much.

Support mkl 2023

Is your feature request related to a problem? Please describe.
The latest version of pyinterp on conda-forge is limited to mkl <2023. However, mkl 2023 is now available.

Describe the solution you'd like
A package available on conda-forge which supports mkl 2023.

Interpolator fails on singleton dimension

Describe the bug

For a 3D variable (time, lat, lon), the interpolation works fine if len(time) > 1, but fails if len(time) == 1 (singleton dimension).

To Reproduce

Code and input/output files bellow. The output file ssh1t_interp.nc only has NaNs.

import numpy as np
import xarray as xr

import pyinterp
import pyinterp.fill
import pyinterp.backends.xarray

# infile = '/tmp/ssh2t.nc'   # OK
infile = '/tmp/ssh1t.nc'   # FAIL

outfile = infile.replace('.nc', '_interp.nc')

ds = xr.open_dataset(infile)

# output coords
time = ds['time']

dx = dy = 1 / 10
latitude = np.arange(ds['latitude'].min(), ds['latitude'].max(), dy)
longitude = np.arange(ds['longitude'].min(), ds['longitude'].max(), dx)

coords = dict(time=time, latitude=latitude, longitude=longitude)

# create output dataset with the new coords
ds_new = xr.Dataset(coords)

mx, my, mt = np.meshgrid(*[ds_new[x].values
                           for x in ['longitude', 'latitude', 'time']],
                         indexing='ij')

for var in ds.data_vars:
    print(f'Interpolating {var}')

    interpolator = pyinterp.backends.xarray.Grid3D(ds[var])
    bicubic = interpolator.bicubic(dict(longitude=mx.flatten(),
                                        latitude=my.flatten(),
                                        time=mt.flatten()))
    bicubic = bicubic.reshape(mx.shape).transpose(2, 1, 0)

    ds_new[var] = (list(ds.coords), bicubic, ds.attrs, ds.encoding)

ds_new.to_netcdf(outfile)
ds.close()

Input/output files:
ssh.zip

Expected behavior

Interpolation for dataarrays with singleton dimension (including time and depth/height dimension) for 3D (time, lat, lon) and 4D (time,depth/height,lat,lon) grids.

Thank you.

openblas package in conda is wrongly linking to mkl

hi,

please fix the conda-forge build:
pyinterp 2023.2.1 openblas_py38hf1160e9_101

it is actually mkl version of the build.
please provide openblas build package to conda-forge.

thanks

`No such file or directory` Error on Pip Install

Describe the bug
Receiving an error when trying to pip install pyinterp

error: [Errno 2] No such file or directory: 'C:\Users\User\AppData\Local\Temp\pip-install-0uuc34iv\pyinterp\third_party\unqlite\unqlite.h'

To Reproduce
pip install pyinterp

**Pyinterp/Numpy/Python version information**

numpy version: 1.19.4
sys version: 3.7.8

Support intel_openmp package

Is your feature request related to a problem? Please describe.
On conda-forge, this package requires the llvm-openmp package, which conflicts with the intel-openmp package.

Describe the solution you'd like
I'd like a version of this package which supports the intel-openmp package, which is required for some packages, such as pytorch.

pyinterp.fill error and error in doc

Errors:

Exception has occurred: AttributeError
module 'pyinterp' has no attribute 'tests'
  File "/root/lastfocus/goes_abi_geo_read.py", line 19, in <module>
    ds = pyi.tests.load_grid2d()
         ^^^^^^^^^
AttributeError: module 'pyinterp' has no attribute 'tests'

To Reproduce
Copy code from Documentation

import pyinterp as pyi
import pyinterp.fill as pf

ds = pyinterp.tests.load_grid2d()
grid = pyinterp.backends.xarray.Grid2D(ds.mss)

Trying to take a 2D data on a regular grid (regional map) that contains some missing values and I would like to use pyinterp.fill.loess to do it. Not sure how the data shape is suppose to be.

large array interpolation, distribution

Let's a consider a situation where the input grid for an interpolation is too large to fit in memory.
A natural way to distribute the interpolation would be to:

1. load chunks of the grid on a distributed cluster (dask of course)
2. perform interpolations for each of the chunks and retain information about distances between target points and input grid.
3. collect results and choose which ones corresponds to the best anwser.

My interrogation is with 2.: it is possible to retain distance information between target points and input grid points that are used for the interpolation?

conda install not working

I cannot intall the library with:

(equinox) [pontea@visu03 drifters]$ conda install pangeo-pyinterp  -c fbriol
Solving environment: failed

PackagesNotFoundError: The following packages are not available from current channels:

  - pangeo-pyinterp

Current channels:

There seems to be a pyinterp library however:

(equinox) [pontea@visu03 drifters]$ conda search -c fbriol | grep pyinterp
pyinterp                  0.0.2  py37ha8d69ae_0  fbriol              

It this what should be used?