Repository: sunpy/ndcube

Issue:
Extend reproject method to support other algorithms
## Outline

The `adaptive` and `exact` reproject algorithms are applicable to image axes only.  This is why PR #434 does not support them.  This issue summarizes the way forward for extending `NDCube.reproject_to` to those other algorithms.
* Initially, we will only support these algorithms in the case of a 2-D celestial `NDCube`.  (Supporting more general cases are discussed further below.)
* If user requests `adaptive` or `exact` to be used then check the original WCS.
  * If not 2-D, raise error
  * If 2-D but not celestial, raise error.
  * Otherwise proceed.

The question now is how to check for celestial axes.

## Checking For Celestial Axes
In order to fast-track this feature, we will use a shortcut for identifying celestial axes which we hope will work for most cases. (See below for the future more generalized approach.)

Check the physical types of the axes against a finite list options.  These will include that the physical type will include:
* `pos.` and `.lat` or `.lon`
* `.ra`
* `.dec`
* Add other cases here.

Because we hope to deal with this more generally in the future, we should check against the UCD1+ vocab list for other valid possibilities.  But beyond that we will encourage users to inform us of other cases which we can then add to our list.  Before this list gets too long we hope to have a more generalized implementation which will remove the need for this checklist.

## Future Generalized Implementation
In general, these algorithms should be valid for any >=2-D cube which contains two dependent spatial axes.  The reprojection could then be done in these two dimensions while keeping the rest unchanged.  To do this we need a way of identifying invariant axes.  This is discussed in #432.  Once that is implemented, we can leverage it remove the above limitations.

## Justification for Immediate Implementation
Supporting other algorithms for 2-D celestial cubes is needed by SunPy Map.  By using the above quick-and-dirty implementation, reproject can be released as part of 2.0 and enable Map to use reproject's full power.  The more generalized subsequent implementation will have no impact on Map's use case.



Discussion:
pep8speaks: Hello @adwaitbhope! Thanks for opening this PR. 

* In the file [`ndcube/ndcube.py`](https://github.com/sunpy/ndcube/blob/bcc423a94958406abe3f32d78e2a78742f810e0a/ndcube/ndcube.py):

> [Line 705:30](https://github.com/sunpy/ndcube/blob/bcc423a94958406abe3f32d78e2a78742f810e0a/ndcube/ndcube.py#L705): [E222](https://duckduckgo.com/?q=pep8%20E222) multiple spaces after operator
> [Line 657:38](https://github.com/sunpy/ndcube/blob/bcc423a94958406abe3f32d78e2a78742f810e0a/ndcube/ndcube.py#L657): [E741](https://duckduckgo.com/?q=pep8%20E741) ambiguous variable name 'l'







Do see the [Hitchhiker's guide to code style](https://goo.gl/hqbW4r)

Patch:
diff --git a/changelog/448.feature.rst b/changelog/448.feature.rst
new file mode 100644
index 0000000..0de881d
--- /dev/null
+++ b/changelog/448.feature.rst
@@ -0,0 +1,1 @@
+Extends `~.NDCube.reproject_to` functionality by supporting ``adaptive`` and ``exact`` algorithms for an `~NDCube` with 2D celestial WCS.
diff --git a/ndcube/ndcube.py b/ndcube/ndcube.py
index 3dc1f97..189fb93 100644
--- a/ndcube/ndcube.py
+++ b/ndcube/ndcube.py
@@ -705,12 +705,33 @@ class NDCubeBase(NDCubeSlicingMixin, NDCubeABC):
         # Creating a new NDCubeSequence with the result_cubes and common axis as axis
         return NDCubeSequence(result_cubes, meta=self.meta)
 
-    def reproject_to(self, target_wcs, shape_out=None, order='bilinear', output_array=None, return_footprint=False):
+    def _validate_algorithm_and_order(self, algorithm, order):
+        order_compatibility = {
+            'interpolation': ['nearest-neighbor', 'bilinear', 'biquadratic', 'bicubic'],
+            'adaptive': ['nearest-neighbor', 'bilinear'],
+            'exact': []
+        }
+
+        if algorithm in order_compatibility:
+            if order_compatibility[algorithm] and order not in order_compatibility[algorithm]:
+                raise ValueError(f"For '{algorithm}' algorithm, the 'order' argument must be "
+                                 f"one of {', '.join(order_compatibility[algorithm])}.")
+
+        else:
+            raise ValueError(f"The 'algorithm' argument must be one of "
+                             f"{', '.join(order_compatibility.keys())}.")
+
+    def reproject_to(self, target_wcs, algorithm='interpolation', shape_out=None, order='bilinear',
+                     output_array=None, parallel=False, return_footprint=False):
         """
         Reprojects this NDCube to the coordinates described by another WCS object.
 
         Parameters
         ----------
+        algorithm: `str`
+            The algorithm to use for reprojecting. This can be any of: 'interpolation', 'adaptive',
+            and 'exact'.
+
         target_wcs : `astropy.wcs.wcsapi.BaseHighLevelWCS`, `astropy.wcs.wcsapi.BaseLowLevelWCS`,
             or `astropy.io.fits.Header`
             The WCS object to which the ``NDCube`` is to be reprojected.
@@ -722,13 +743,24 @@ class NDCubeBase(NDCubeSlicingMixin, NDCubeABC):
             (if available) from the low level API of the ``target_wcs`` is used.
 
         order: `int` or `str`
-            The order of the interpolation. This can be any of: 'nearest-neighbor', 'bilinear',
-            'biquadratic', or 'bicubic'.
+            The order of the interpolation (used only when the 'interpolation' or 'adaptive'
+            algorithm is selected).
+            For 'interpolation' algorithm, this can be any of: 'nearest-neighbor', 'bilinear',
+            'biquadratic', and 'bicubic'.
+            For 'adaptive' algorithm, this can be either 'nearest-neighbor' or 'bilinear'.
 
         output_array: `numpy.ndarray`, optional
             An array in which to store the reprojected data. This can be any numpy array
             including a memory map, which may be helpful when dealing with extremely large files.
 
+        parallel: `bool` or `int`
+            Flag for parallel implementation (used only when the 'exact' algorithm is selected).
+            If ``True``, a parallel implementation is chosen and the number of processes is
+            selected automatically as the number of logical CPUs detected on the machine.
+            If ``False``, a serial implementation is chosen.
+            If the flag is a positive integer n greater than one, a parallel implementation
+            using n processes is chosen.
+
         return_footprint: `bool`
             Whether to return the footprint in addition to the output NDCube.
 
@@ -747,18 +779,28 @@ class NDCubeBase(NDCubeSlicingMixin, NDCubeABC):
         However, ``meta`` and ``global_coords`` are copied to the output ``NDCube``.
         """
         try:
-            from reproject import reproject_interp
+            from reproject import reproject_adaptive, reproject_exact, reproject_interp
+            from reproject.wcs_utils import has_celestial
         except ModuleNotFoundError:
             raise ImportError("The NDCube.reproject_to method requires the optional package `reproject`.")
 
         if isinstance(target_wcs, Mapping):
             target_wcs = WCS(header=target_wcs)
 
+        low_level_target_wcs = utils.wcs.get_low_level_wcs(target_wcs, 'target_wcs')
+
+        # 'adaptive' and 'exact' algorithms work only on 2D celestial WCS.
+        if algorithm == 'adaptive' or algorithm == 'exact':
+            if low_level_target_wcs.pixel_n_dim != 2 or low_level_target_wcs.world_n_dim != 2:
+                raise ValueError('For adaptive and exact algorithms, target_wcs must be 2D.')
+
+            if not has_celestial(target_wcs):
+                raise ValueError('For adaptive and exact algorithms, '
+                                 'target_wcs must contain celestial axes only.')
+
         if not utils.wcs.compare_wcs_physical_types(self.wcs, target_wcs):
             raise ValueError('Given target_wcs is not compatible with this NDCube, the physical types do not match.')
 
-        low_level_target_wcs = utils.wcs.get_low_level_wcs(target_wcs, 'target_wcs')
-
         # If shape_out is not specified explicity, try to extract it from the low level WCS
         if not shape_out:
             if hasattr(low_level_target_wcs, 'array_shape') and low_level_target_wcs.array_shape is not None:
@@ -766,10 +808,20 @@ class NDCubeBase(NDCubeSlicingMixin, NDCubeABC):
             else:
                 raise ValueError("shape_out must be specified if target_wcs does not have the array_shape attribute.")
 
-        data = reproject_interp(self, output_projection=target_wcs,
-                                shape_out=shape_out, order=order,
-                                output_array=output_array,
-                                return_footprint=return_footprint)
+        self._validate_algorithm_and_order(algorithm, order)
+
+        if algorithm == 'interpolation':
+            data = reproject_interp(self, output_projection=target_wcs, shape_out=shape_out,
+                                    order=order, output_array=output_array,
+                                    return_footprint=return_footprint)
+
+        elif algorithm == 'adaptive':
+            data = reproject_adaptive(self, output_projection=target_wcs, shape_out=shape_out,
+                                      order=order, return_footprint=return_footprint)
+
+        elif algorithm == 'exact':
+            data = reproject_exact(self, output_projection=target_wcs, shape_out=shape_out,
+                                   parallel=parallel, return_footprint=return_footprint)
 
         if return_footprint:
             data, footprint = data


Test patch:
diff --git a/ndcube/tests/test_ndcube.py b/ndcube/tests/test_ndcube.py
index 708fbc8..1666b91 100644
--- a/ndcube/tests/test_ndcube.py
+++ b/ndcube/tests/test_ndcube.py
@@ -566,13 +566,13 @@ def test_initialize_from_ndcube(ndcube_3d_l_ln_lt_ectime):
     assert ec is not ec3
 
 
-def test_reproject(ndcube_4d_ln_l_t_lt, wcs_4d_lt_t_l_ln):
+def test_reproject_interpolation(ndcube_4d_ln_l_t_lt, wcs_4d_lt_t_l_ln):
     target_wcs_header = wcs_4d_lt_t_l_ln.low_level_wcs.to_header()
     target_wcs_header['CDELT3'] = 0.1   # original value = 0.2
     target_wcs = astropy.wcs.WCS(header=target_wcs_header)
     shape_out = (5, 20, 12, 8)
 
-    resampled_cube = ndcube_4d_ln_l_t_lt.reproject_to(target_wcs, shape_out)
+    resampled_cube = ndcube_4d_ln_l_t_lt.reproject_to(target_wcs, shape_out=shape_out)
 
     assert ndcube_4d_ln_l_t_lt.data.shape == (5, 10, 12, 8)
     assert resampled_cube.data.shape == (5, 20, 12, 8)
@@ -582,14 +582,14 @@ def test_reproject_invalid_wcs(ndcube_4d_ln_l_t_lt, wcs_3d_lt_ln_l):
     shape_out = (5, 20, 12, 8)
 
     with pytest.raises(Exception):
-        _ = ndcube_4d_ln_l_t_lt.reproject_to(wcs_3d_lt_ln_l, shape_out)
+        _ = ndcube_4d_ln_l_t_lt.reproject_to(wcs_3d_lt_ln_l, shape_out=shape_out)
 
 
 def test_reproject_with_header(ndcube_4d_ln_l_t_lt, wcs_4d_lt_t_l_ln):
     target_wcs_header = wcs_4d_lt_t_l_ln.low_level_wcs.to_header()
     shape_out = (5, 20, 12, 8)
 
-    _ = ndcube_4d_ln_l_t_lt.reproject_to(target_wcs_header, shape_out)
+    _ = ndcube_4d_ln_l_t_lt.reproject_to(target_wcs_header, shape_out=shape_out)
 
 
 def test_reproject_return_footprint(ndcube_4d_ln_l_t_lt, wcs_4d_lt_t_l_ln):
@@ -598,7 +598,7 @@ def test_reproject_return_footprint(ndcube_4d_ln_l_t_lt, wcs_4d_lt_t_l_ln):
     target_wcs = astropy.wcs.WCS(header=target_wcs_header)
     shape_out = (5, 20, 12, 8)
 
-    resampled_cube, footprint = ndcube_4d_ln_l_t_lt.reproject_to(target_wcs, shape_out,
+    resampled_cube, footprint = ndcube_4d_ln_l_t_lt.reproject_to(target_wcs, shape_out=shape_out,
                                                                  return_footprint=True)
 
     assert ndcube_4d_ln_l_t_lt.data.shape == (5, 10, 12, 8)
@@ -630,3 +630,72 @@ def test_wcs_type_after_init(ndcube_3d_ln_lt_l, wcs_3d_l_lt_ln):
     cube = NDCube(ndcube_3d_ln_lt_l.data[slices], low_level_wcs)
     # Check the WCS has been converted to high level but NDCube init.
     assert isinstance(cube.wcs, BaseHighLevelWCS)
+
+
+def test_reproject_adaptive(ndcube_2d_ln_lt, wcs_2d_lt_ln):
+    shape_out = (10, 12)
+    resampled_cube = ndcube_2d_ln_lt.reproject_to(wcs_2d_lt_ln, algorithm='adaptive',
+                                                  shape_out=shape_out)
+
+    assert ndcube_2d_ln_lt.data.shape == (10, 12)
+    assert resampled_cube.data.shape == (10, 12)
+
+
+def test_reproject_exact(ndcube_2d_ln_lt, wcs_2d_lt_ln):
+    shape_out = (10, 12)
+    resampled_cube = ndcube_2d_ln_lt.reproject_to(wcs_2d_lt_ln, algorithm='exact',
+                                                  shape_out=shape_out)
+
+    assert ndcube_2d_ln_lt.data.shape == (10, 12)
+    assert resampled_cube.data.shape == (10, 12)
+
+
+def test_reproject_invalid_algorithm(ndcube_4d_ln_l_t_lt, wcs_4d_lt_t_l_ln):
+    with pytest.raises(ValueError):
+        _ = ndcube_4d_ln_l_t_lt.reproject_to(wcs_4d_lt_t_l_ln, algorithm='my_algorithm',
+                                             shape_out=(5, 10, 12, 8))
+
+
+def test_reproject_invalid_order(ndcube_2d_ln_lt, wcs_2d_lt_ln):
+    shape_out = (10, 12)
+
+    # Supported for interpolation
+    for order in ['nearest-neighbor', 'bilinear', 'biquadratic', 'bicubic']:
+        _ = ndcube_2d_ln_lt.reproject_to(wcs_2d_lt_ln, algorithm='interpolation',
+                                         shape_out=shape_out, order=order)
+
+    # Supported for adaptive
+    for order in ['nearest-neighbor', 'bilinear']:
+        _ = ndcube_2d_ln_lt.reproject_to(wcs_2d_lt_ln, algorithm='adaptive',
+                                         shape_out=shape_out, order=order)
+
+    # Not supported for adaptive
+    for order in ['biquadratic', 'bicubic', 'my_order']:
+        with pytest.raises(ValueError):
+            _ = ndcube_2d_ln_lt.reproject_to(wcs_2d_lt_ln, algorithm='adaptive',
+                                             shape_out=shape_out, order=order)
+
+    # 'exact' does not need 'order'
+    _ = ndcube_2d_ln_lt.reproject_to(wcs_2d_lt_ln, algorithm='exact', shape_out=shape_out)
+
+
+def test_reproject_adaptive_incompatible_wcs(ndcube_4d_ln_l_t_lt, wcs_4d_lt_t_l_ln,
+                                             wcs_1d_l, ndcube_1d_l):
+    with pytest.raises(ValueError):
+        _ = ndcube_1d_l.reproject_to(wcs_1d_l, algorithm='adaptive',
+                                     shape_out=(10,))
+
+    with pytest.raises(ValueError):
+        _ = ndcube_4d_ln_l_t_lt.reproject_to(wcs_4d_lt_t_l_ln, algorithm='adaptive',
+                                             shape_out=(5, 10, 12, 8))
+
+
+def test_reproject_exact_incompatible_wcs(ndcube_4d_ln_l_t_lt, wcs_4d_lt_t_l_ln,
+                                          wcs_1d_l, ndcube_1d_l):
+    with pytest.raises(ValueError):
+        _ = ndcube_1d_l.reproject_to(wcs_1d_l, algorithm='exact',
+                                     shape_out=(10,))
+
+    with pytest.raises(ValueError):
+        _ = ndcube_4d_ln_l_t_lt.reproject_to(wcs_4d_lt_t_l_ln, algorithm='exact',
+                                             shape_out=(5, 10, 12, 8))