Compute residuals to the timeseries inversion, save to raster outputs (…

…#523)

* Add `censored_lstsq` to solve least squares with missing data

* Put conncomp file list first

* [skip ci] remove dup line

* Fix `invert_unw_stack` to set output pixel `nodata` based on unw nodata

* start writing residuals to a raster

* Fix bug in gdal overview threadd lzw creation

When doing `GDAL_NUM_THREADS=2` and LZW, this came up
  File "/Users/staniewi/repos/dolphin/src/dolphin/_overviews.py", line 95, in create_image_overviews
    ds.BuildOverviews(resampling.value, levels)
  File "/Users/staniewi/miniforge3/envs/mapping-312/lib/python3.12/site-packages/osgeo/gdal.py", line 4378, in BuildOverviews
    return _gdal.Dataset_BuildOverviews(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: /vsimem/decompress_0x312935528.tif:Using code not yet in table

After turning the python multithreading off, and the GDAL_NUM_THREADS=1,
it was still giving this:

RuntimeError: 20210606_20210618.tif, band 1: IReadBlock failed at X offset 0, Y offset 2: TIFFReadEncodedTile() failed.
May be caused by: TIFFReadEncodedTile() failed.
May be caused by: /Users/staniewi/repos/dolphin/docs/notebooks/work-walkthrough/timeseries/20210606_20210618.tif:Using code not yet in table

This was on a Mac with Gdal 3.9.0
$ gdalinfo --version
GDAL 3.9.3, released 2024/10/07

* Add a `set_raster_nodata` to match the `get_`

* add missing io tests

* Fix occasional sorting issue for residuals, nodata outputs

* remove extra utils import

* Get residuals per date working

* fix overview by setting nodata value in `GA_Update` mode

* setup residual path tests, add to workflows and `extra_reference` logic

* fix tests

* fix last resid path problem

src/dolphin/_overviews.py

@@ -54,6 +54,7 @@ def create_image_overviews(
     resampling: Resampling | None = None,
     external: bool = False,
     compression: str = "LZW",
+    num_gdal_threads: int = 1,
 ):
     """Add GDAL compressed overviews to an existing file.
 
@@ -72,9 +73,13 @@ def create_image_overviews(
         if not specifying `image_type`.
     external : bool, default = False
         Use external overviews (.ovr files).
-    compression: str, default = "LZW"
+    compression: str, default = "DEFLATE"
         Compression algorithm to use for overviews.
         See https://gdal.org/programs/gdaladdo.html for options.
+    num_gdal_threads: int, default = 1
+        Number of threads to use for overview creation.
+        For more details, see https://gdal.org/programs/gdaladdo.html and
+        https://gdal.org/en/stable/user/configoptions.html#config-GDAL_NUM_THREADS
 
     """
     if image_type is None and resampling is None:
@@ -91,7 +96,7 @@ def create_image_overviews(
         return
 
     gdal.SetConfigOption("COMPRESS_OVERVIEW", compression)
-    gdal.SetConfigOption("GDAL_NUM_THREADS", "2")
+    gdal.SetConfigOption("GDAL_NUM_THREADS", str(num_gdal_threads))
     ds.BuildOverviews(resampling.value, levels)
 
 

src/dolphin/io/_core.py

@@ -50,6 +50,7 @@
     "load_gdal",
     "set_raster_description",
     "set_raster_metadata",
+    "set_raster_nodata",
     "set_raster_units",
     "write_arr",
     "write_block",
@@ -83,9 +84,10 @@
 }
 
 
-def _get_gdal_ds(filename: Filename) -> gdal.Dataset:
+def _get_gdal_ds(filename: Filename, update: bool = False) -> gdal.Dataset:
+    mode = gdal.GA_Update if update else gdal.GA_ReadOnly
     if str(filename).startswith("s3://"):
-        return gdal.Open(S3Path(str(filename)).to_gdal())
+        return gdal.Open(S3Path(str(filename)).to_gdal(), mode)
     return gdal.Open(fspath(filename))
 
 
@@ -284,6 +286,29 @@ def get_raster_nodata(filename: Filename, band: int = 1) -> Optional[float]:
     return ds.GetRasterBand(band).GetNoDataValue()
 
 
+def set_raster_nodata(filename: Filename, nodata: float, band: int | None = None):
+    """Set the nodata value for a raster.
+
+    Parameters
+    ----------
+    filename : Filename
+        Path to the file to load.
+    nodata : float
+        The nodata value to set.
+    band : int, optional
+        The band to set the nodata value for, by default None
+        (sets the nodata value for all bands).
+
+    """
+    ds = gdal.Open(fspath(filename), gdal.GA_Update)
+    if band is None:
+        for i in range(ds.RasterCount):
+            ds.GetRasterBand(i + 1).SetNoDataValue(nodata)
+    else:
+        ds.GetRasterBand(band).SetNoDataValue(nodata)
+    ds = None
+
+
 def get_raster_crs(filename: Filename) -> CRS:
     """Get the CRS from a file.