[WIP] Multiprocessing interpolation and reprojection

[vschaffn]

Resolves #648.

Description

This PR introduces raster interpolation and reprojection using multiprocessing to optimize memory usage and improve performance.

Code changes

• Multiprocessing interpolation : The function multiproc_interp_points leverages multiprocessing to create a separate process for each point to be interpolated. For each point, the function calculates the smallest possible raster window to open based on the interpolation method and point coordinates, and use the raster.crop method to open the window. The raster.interpolate method is then applied to the cropped raster window, allowing for efficient interpolation without loading the full raster into memory.
• Multiprocessing reprojection: The function multiproc_reproject utilizes raster tiling (see Implement Tiling for Multiprocessing #652) to perform multiprocessing-based reprojection. A separate process is created for each tile of the raster. Each process opens the corresponding tile using the raster.icrop method, calculates the bounds of the reprojected tile, snaps the bounds to the target reprojected grid, and reprojects the tile using the raster.reproject method. The reprojected tiles are written separately to disk, with safeguards in place to prevent data overwriting.

Tests

• Unit tests have been implemented to verify the correctness of both multiproc_interp_points and multiproc_reproject. The test results are compared against the behavior of the raster.interpolate and raster.reproject methods to ensure consistency.

Note:

Currently, there are some tile alignment issues when performing reprojection operations that involve more than just translation. Further investigation is required to address these challenges.

Difference between tiled reprojection and classic reprojection (exploradores_aster_dem):

• Changed CRS:
  

[rhugonnet]

Hi @vschaffn @adebardo, nice to see the progress!

Summary: Upon reading the code in delayed_multiproc, I'm fairly sure that in many cases the reproject won't work, and that interp might work but won't be computationally optimal.

Why?
For reproject, it's due to shape deformation during CRS reprojection, you can't work simply on extents. You need the equivalent of the GeoGrid classes, a vectorized format of the grid shape. Because of the nature of the problem, writing these steps as embedded functions would probably take three times the code size including many loops and be hard to grasp, which is why I wrote it this way (and they also did in other packages, like OpenDataCube-Geo).
For interp, it's because loading many tiny windows will be extremely slow. It's better to load a whole block that fits in memory, and perform interpolation on all the points of this block at once.

Moving forward:
Joining my previous comment (#655 (comment)), I suggest to simply re-use all the internal steps of dask_reproject and dask_interp functions that were already developed and extensively tested, and only replace the bits about loading/tiling that work differently without Dask and using Multiprocessing instead.
For example, in dask_reproject, nothing uses Dask until this line:

geoutils/geoutils/raster/delayed.py

Line 783 in b4da49b

blocks = darr.to_delayed().ravel()

So all the projected grid matching calculations to map the input chunks to output chunks written above that line (using the GeoGrid classes I mentioned) can simply stay the same 😉

[rhugonnet]

Also note that Rasterio has reprojection errors that can appear (especially when comparing a tiled reprojection to a reprojection on the full raster), due to internal inconsistencies...
But those look stripped, like this: rasterio/rasterio#2052 (comment)

[rhugonnet]

To clarify one more remark: In multiprocessing_reproject(), you'd have to slightly change the logic you currently have within your tile loop (line 404).

Due to shape deformations, the mapping of input-output tiles between any CRS requires more than an overlap value, and also depends on the input/output chunksizes (what the GeoGrid code mentioned above does).
So you'll have to loop for each output tile (always 1 by 1), and every time open a combination of input tiles (could be 1, 2, 3, 4, 5, ...).

You can essentially turn this list comprehension here into your loop for multiprocessing_reproject:

geoutils/geoutils/raster/delayed.py

Line 784 in b4da49b

# Run the delayed reprojection, looping for each destination block

And let reproject_block stick the pieces of input tiles in the right places (you might not have a perfect square with all the tiles you opened), and give you the output:

geoutils/geoutils/raster/delayed.py

Line 644 in b4da49b

def _delayed_reproject_per_block(

[rhugonnet]

I should have thought about this last week (I was mostly off work and only popped for the comment above, so I didn't think about the big picture, sorry!):
As I just mentioned in GlacioHack/xdem#699, the current code for reproject developed by @vschaffn will still be extremely useful, as it is essentially the multiprocessing equivalent of dask.map_overlap that we need absolutely everywhere (all terrain attributes such as the slope, all kernel-based filters, re-gridding in Coreg.apply(), and much more 😊).

It is only reproject that is a tricky exception! 😅