Merge pull request #25 from FrontierDevelopmentLab/bugfix/corrupted-bands

Fix rescaling bug

Author: frandorr

.gcloudignore

@@ -4,7 +4,8 @@ tests
 *.pyo
 *.pyd
 notebooks
-outputs
+output
+images
 docs
 htmlcov
 dist
@@ -14,6 +15,8 @@ ci
 .vscode
 .mypy_cache
 .ipynb_checkpoints
-uis
 __pycache__
 .pytest_cache
+.hydra
+ui
+*.ipynb

providers/gcp/main.py

@@ -106,7 +106,6 @@ def extract_patches():
         # common data
         storage_gs_path = request_json["storage_gs_path"]
         bands = request_json["bands"]
-        resolution = request_json["resolution"]
         job_id = request_json["job_id"]
 
         fs = gcsfs.GCSFileSystem()
@@ -150,16 +149,16 @@ def extract_patches():
                 "Environment variable MONITOR_TABLE not set. Unable to push task status to Monitor",
             )
 
+        archive_resolution = int(
+            min([b["gsd"] for _, b in BAND_INFO[constellation].items()]),
+        )
+
         patches = task_mosaic_patches(
             cloud_fs=fs,
             download_f=download_blob,
             task=task,
             method="first",
-            resolution=resolution,
-        )
-
-        archive_resolution = int(
-            min([b["gsd"] for kk, b in BAND_INFO[constellation].items()]),
+            resolution=archive_resolution,
         )
 
         logger.info(f"Ready to store {len(patches)} patches at {storage_gs_path}.")
@@ -169,7 +168,6 @@ def extract_patches():
             patches,
             task,
             bands,
-            resolution,
             archive_resolution,
         )
 

src/satextractor/builder/gcp_builder.py

@@ -150,7 +150,7 @@ def build_docker_image(self):
             "--tag",
             f"{self.image_region_code}/{self.project}/{self.user_id}-stacextractor",
         ]
-
+        logger.info(cmd)
         p = run(cmd, text=True, stdout=subprocess.DEVNULL)
         p.check_returncode()
 

src/satextractor/deployer/gcp_deployer.py

@@ -29,7 +29,6 @@ def deploy_tasks(credentials, extraction_tasks, storage_path, chunk_size, topic)
             extraction_task=extraction_task_data,
             bands=list(BAND_INFO[task.constellation].keys()),
             chunks=(1, 1, chunk_size, chunk_size),
-            resolution=int(BAND_INFO[task.constellation][task.band]["gsd"]),
         )
         data = json.dumps(data, default=str)
 

src/satextractor/extractor/extractor.py

@@ -12,6 +12,7 @@
 from osgeo import osr
 from rasterio import warp
 from rasterio.crs import CRS
+from rasterio.enums import Resampling
 from rasterio.merge import merge as riomerge
 from satextractor.models import ExtractionTask
 from satextractor.models import Tile
@@ -106,6 +107,7 @@ def download_and_extract_tiles_window_COG(
     # set the transforms for the output file
     dst_transform = Affine(resolution, 0.0, left, 0.0, -resolution, top)
     out_shp = (int((right - left) / resolution), int((top - bottom) / resolution))
+    logger.debug(f"Affine with resolution: {dst_transform} and out_shp {out_shp} ")
 
     outfiles = []
 
@@ -136,6 +138,7 @@ def download_and_extract_tiles_window_COG(
             transform=dst_transform,
             crs=CRS.from_epsg(epsg),
             dtype=rst_arr.dtype,
+            resampling=Resampling.bilinear,
         ) as dst:
 
             dst.write(rst_arr, indexes=1)
@@ -196,8 +199,8 @@ def download_and_extract_tiles_window(
             projWin=proj_win,
             projWinSRS=f"EPSG:{epsg}",
             xRes=resolution,
-            yRes=resolution,
-            resampleAlg="cubic",
+            yRes=-resolution,
+            resampleAlg="bilinear",
             creationOptions=["QUALITY=100", "REVERSIBLE=YES"],
         )
         file = None

src/satextractor/storer/storer.py

@@ -5,7 +5,6 @@
 import numpy as np
 import zarr
 from satextractor.models import ExtractionTask
-from scipy.ndimage import zoom
 
 
 def store_patches(
@@ -14,7 +13,6 @@ def store_patches(
     patches: List[np.ndarray],
     task: ExtractionTask,
     bands: List[str],
-    patch_resolution: int,
     archive_resolution: int,
 ):
     """Store a list of patches in storage path.
@@ -48,11 +46,6 @@ def store_patches(
         timestamp_idx = timestamps.index(task.sensing_time)
         patch = patches[i]
 
-        # maybe resize -> bicubic upsample
-        if patch_resolution != archive_resolution:
-            patch = zoom(patch, int(patch_resolution / archive_resolution), order=3)
-
-        # in patch resolution
         if patch.shape != size:
             pad_x = int(size[0] - patch.shape[0])
             pad_y = int(size[1] - patch.shape[1])