Repository: ASFHyP3/hyp3-autorift

Issue:
Landsat 7 + 8 pairs may not be handled well
Since we only look at the reference scene to [determine the platform](https://github.com/ASFHyP3/hyp3-autorift/blob/develop/hyp3_autorift/process.py#L319), we may have some issues with the secondary scene:
- [x] L7 reference scene w/ a L8 secondary will attempt to apply the Wallis filter and bonk:
    https://github.com/ASFHyP3/hyp3-autorift/blob/develop/hyp3_autorift/process.py#L269-L282
    
- [x] L8 reference scene w/ a L7 secondary will not have the Wallis filter applied to it:
    https://github.com/ASFHyP3/hyp3-autorift/blob/develop/hyp3_autorift/process.py#L368-L370

Fortunately, in our current campaign pair list, 4/5 pairs are only crossed with 4/5 pairs, and 9 pairs are only crossed with 8, so we should only see issues around 7+8 pairs.

```python
>>> import geopandas as gpd
>>> df = gpd.read_parquet('l45789.parquet')
>>> df.groupby(['ref_mission', 'sec_mission']).reference.count()
ref_mission  sec_mission
L4           L4                1058
             L5                3203
L5           L4                2988
             L5              767024
L7           L7             2003704
             L8              474339
L8           L7              416539
             L8             4103107
             L9              367178
L9           L8              112403
             L9              109936
```

**Fix:** AutoRIFT pre-processing defaults to the strictest filtering, so we should allow L8 scenes to be wallis filled if either pair is L7. And we should check both ref, sec platform when deciding to run the pre-processing filters or not in `process.py`

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index e092e35..a9108fb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,11 @@ and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/)
 and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 
+## [0.10.4]
+
+### Fixed
+* Landsat 7+8 pairs will be filtered appropriately; see [#201](https://github.com/ASFHyP3/hyp3-autorift/issues/201)
+
 ## [0.10.3]
 
 ### Added
diff --git a/hyp3_autorift/io.py b/hyp3_autorift/io.py
index d7a1b22..9876b6d 100644
--- a/hyp3_autorift/io.py
+++ b/hyp3_autorift/io.py
@@ -168,7 +168,7 @@ def load_geospatial(infile: str, band: int = 1):
 
 
 def write_geospatial(outfile: str, data, transform, projection, nodata,
-                     driver: str = 'GTiff', dtype: int = gdal.GDT_Float64):
+                     driver: str = 'GTiff', dtype: int = gdal.GDT_Float64) -> str:
     driver = gdal.GetDriverByName(driver)
 
     rows, cols = data.shape
diff --git a/hyp3_autorift/process.py b/hyp3_autorift/process.py
index a56b7f6..b887d5e 100644
--- a/hyp3_autorift/process.py
+++ b/hyp3_autorift/process.py
@@ -11,7 +11,7 @@ import xml.etree.ElementTree as ET
 from datetime import datetime
 from pathlib import Path
 from secrets import token_hex
-from typing import Optional, Tuple
+from typing import Callable, Optional, Tuple
 
 import boto3
 import botocore.exceptions
@@ -227,12 +227,11 @@ def get_s1_primary_polarization(granule_name):
     raise ValueError(f'Cannot determine co-polarization of granule {granule_name}')
 
 
-def create_filtered_filepath(path: str):
+def create_filtered_filepath(path: str) -> Path:
     parent = (Path.cwd() / 'filtered').resolve()
     parent.mkdir(exist_ok=True)
 
-    out_path = parent / Path(path).name
-    return str(out_path)
+    return parent / Path(path).name
 
 
 def prepare_array_for_filtering(array: np.ndarray, nodata: int) -> Tuple[np.ndarray, np.ndarray]:
@@ -266,32 +265,52 @@ def apply_wallis_nodata_fill_filter(array: np.ndarray, nodata: int) -> Tuple[np.
     return filtered, zero_mask
 
 
-def apply_landsat_filtering(image_path: str, image_platform: str) -> Tuple[Path, Optional[Path]]:
+def _apply_filter_function(image_path: str, filter_function: Callable) -> Tuple[Path, Optional[Path]]:
     image_array, image_transform, image_projection, image_nodata = io.load_geospatial(image_path)
     image_array = image_array.astype(np.float32)
 
+    image_filtered, zero_mask = filter_function(image_array, image_nodata)
+
+    image_new_path = create_filtered_filepath(image_path)
+    _ = io.write_geospatial(str(image_new_path), image_filtered, image_transform, image_projection,
+                            nodata=0, dtype=gdal.GDT_Float32)
+
+    zero_path = None
+    if zero_mask is not None:
+        zero_path = create_filtered_filepath(f'{image_new_path.stem}_zeroMask{image_new_path.suffix}')
+        _ = io.write_geospatial(str(zero_path), zero_mask, image_transform, image_projection,
+                                nodata=np.iinfo(np.uint8).max, dtype=gdal.GDT_Byte)
+
+    return image_new_path, zero_path
+
+
+def apply_landsat_filtering(reference: str, secondary: str) -> Tuple[Path, Optional[Path], Path, Optional[Path]]:
+    reference_platform = get_platform(reference)
+    secondary_platform = get_platform(secondary)
+    if reference_platform > 'L7' and secondary_platform > 'L7':
+        raise NotImplementedError(
+            f'{reference_platform}+{secondary_platform} pairs should be highpass filtered in autoRIFT instead'
+        )
+
     platform_filter_dispatch = {
         'L4': apply_fft_filter,
         'L5': apply_fft_filter,
         'L7': apply_wallis_nodata_fill_filter,
+        'L8': apply_wallis_nodata_fill_filter,  # sometimes paired w/ L7 scenes, so use same filter
     }
-
     try:
-        image_filtered, zero_mask = platform_filter_dispatch[image_platform](image_array, image_nodata)
+        reference_filter = platform_filter_dispatch[reference_platform]
+        secondary_filter = platform_filter_dispatch[secondary_platform]
     except KeyError:
-        raise NotImplementedError(f'Unknown pre-processing filter for satellite platform: {image_platform}')
+        raise NotImplementedError('Unknown pre-processing filter for satellite platform')
 
-    image_new_path = create_filtered_filepath(image_path)
-    image_path = io.write_geospatial(image_new_path, image_filtered, image_transform, image_projection,
-                                     nodata=0, dtype=gdal.GDT_Float32)
+    if reference_filter != secondary_filter:
+        raise NotImplementedError('AutoRIFT not available for image pairs with different preprocessing methods')
 
-    zero_path = None
-    if zero_mask is not None:
-        zero_path = create_filtered_filepath(f'{Path(image_path).stem}_zeroMask{Path(image_path).suffix}')
-        _ = io.write_geospatial(str(zero_path), zero_mask, image_transform, image_projection,
-                                nodata=np.iinfo(np.uint8).max, dtype=gdal.GDT_Byte)
+    reference_path, reference_zero_path = _apply_filter_function(reference, reference_filter)
+    secondary_path, secondary_zero_path = _apply_filter_function(secondary, secondary_filter)
 
-    return image_path, zero_path
+    return reference_path, reference_zero_path, secondary_path, secondary_zero_path
 
 
 def process(reference: str, secondary: str, parameter_file: str = DEFAULT_PARAMETER_FILE,
@@ -365,9 +384,10 @@ def process(reference: str, secondary: str, parameter_file: str = DEFAULT_PARAME
         secondary_metadata = get_lc2_metadata(secondary)
         secondary_path = get_lc2_path(secondary_metadata)
 
-        if platform in ('L4', 'L5', 'L7'):
-            reference_path, reference_zero_path = apply_landsat_filtering(reference_path, platform)
-            secondary_path, secondary_zero_path = apply_landsat_filtering(secondary_path, platform)
+        filter_platform = min([platform, get_platform(secondary)])
+        if filter_platform in ('L4', 'L5', 'L7'):
+            reference_path, reference_zero_path, secondary_path, secondary_zero_path = \
+                apply_landsat_filtering(reference_path, secondary_path)
 
         if reference_metadata['properties']['proj:epsg'] != secondary_metadata['properties']['proj:epsg']:
             log.info('Reference and secondary projections are different! Reprojecting.')
@@ -376,9 +396,6 @@ def process(reference: str, secondary: str, parameter_file: str = DEFAULT_PARAME
             if reference_zero_path and secondary_zero_path:
                 _, _ = io.ensure_same_projection(reference_zero_path, secondary_zero_path)
 
-            elif not isinstance(reference_zero_path, type(secondary_zero_path)):
-                raise NotImplementedError('AutoRIFT not available for image pairs with different preprocessing methods')
-
             reference_path, secondary_path = io.ensure_same_projection(reference_path, secondary_path)
 
         bbox = reference_metadata['bbox']


Test patch:
diff --git a/tests/test_process.py b/tests/test_process.py
index eae954f..42d58d1 100644
--- a/tests/test_process.py
+++ b/tests/test_process.py
@@ -1,5 +1,6 @@
 import io
 from datetime import datetime
+from pathlib import Path
 from re import match
 from unittest import mock
 from unittest.mock import MagicMock, patch
@@ -358,3 +359,67 @@ def test_get_s1_primary_polarization():
         process.get_s1_primary_polarization('S1A_IW_SLC__1SVH_20150706T015744_20150706T015814_006684_008EF7_9B69')
     with pytest.raises(ValueError):
         process.get_s1_primary_polarization('S1A_IW_GRDH_1SVV_20150706T015720_20150706T015749_006684_008EF7_54BA')
+
+
+def test_apply_landsat_filtering(monkeypatch):
+    def mock_apply_filter_function(scene, _):
+        if process.get_platform(scene) < 'L7':
+            return Path(scene), None
+        return Path(scene), Path('zero_mask')
+
+    monkeypatch.setattr(process, '_apply_filter_function', mock_apply_filter_function)
+
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LC09', 'LC09')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LC09', 'LC08')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LC09', 'LE07')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LC09', 'LT05')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LC09', 'LT04')
+
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LC08', 'LC09')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LC08', 'LC08')
+    assert process.apply_landsat_filtering('LC08', 'LE07') == \
+           (Path('LC08'), Path('zero_mask'), Path('LE07'), Path('zero_mask'))
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LC08', 'LT05')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LC08', 'LT04')
+
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LE07', 'LC09')
+    assert process.apply_landsat_filtering('LE07', 'LC08') == \
+           (Path('LE07'), Path('zero_mask'), Path('LC08'), Path('zero_mask'))
+    assert process.apply_landsat_filtering('LE07', 'LE07') == \
+           (Path('LE07'), Path('zero_mask'), Path('LE07'), Path('zero_mask'))
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LE07', 'LT05')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LE07', 'LT04')
+
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LT05', 'LC09')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LT05', 'LC08')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LT05', 'LE07')
+    assert process.apply_landsat_filtering('LT05', 'LT05') == \
+           (Path('LT05'), None, Path('LT05'), None)
+    assert process.apply_landsat_filtering('LT05', 'LT04') == \
+           (Path('LT05'), None, Path('LT04'), None)
+
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LT04', 'LC09')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LT04', 'LC08')
+    with pytest.raises(NotImplementedError):
+        process.apply_landsat_filtering('LT04', 'LE07')
+    assert process.apply_landsat_filtering('LT04', 'LT05') == \
+           (Path('LT04'), None, Path('LT05'), None)
+    assert process.apply_landsat_filtering('LT04', 'LT04') == \
+           (Path('LT04'), None, Path('LT04'), None)