New feature (WIT Charts) in dea-conflux and major package upgrades (for code review and discussion, not ready to merge)

.github/workflows/test.yml

@@ -19,11 +19,6 @@ jobs:
     - name: Pre-pull layers
       run: docker compose pull
 
-    - name: Activate Docker cache
-      uses: satackey/action-docker-layer-caching@v0.0.8
-      # Ignore the failure of a step and avoid terminating the job.
-      continue-on-error: true
-
     - name: Build dea-conflux image
       run: |
         docker compose build

Dockerfile

@@ -1,4 +1,4 @@
-FROM osgeo/gdal:ubuntu-small-3.4.1
+FROM ghcr.io/osgeo/gdal:ubuntu-small-3.6.3
 
 ENV DEBIAN_FRONTEND=noninteractive \
     LC_ALL=C.UTF-8 \
@@ -15,7 +15,7 @@ RUN apt-get update && \
       wget \
       unzip \
       python3-pip \
-      libpq-dev python-dev \
+      libpq-dev python3 \
     && apt-get autoclean && \
     apt-get autoremove && \
     rm -rf /var/lib/{apt,dpkg,cache,log}

constraints.txt

@@ -1,75 +1,3 @@
-affine==2.3.0
-attrs==21.2.0
-boto3==1.20.24
-botocore==1.23.24
-cachetools==4.2.4
-certifi==2021.10.8
-cffi==1.15.1
-cftime==1.5.1.1
-charset-normalizer==3.2.0
-click==8.0.3
-click-plugins==1.1.1
-cligj==0.7.2
-cloudpickle==2.0.0
-coverage==7.3.1
-cryptography==41.0.3
-dask==2021.12.0
-distributed==2021.12.0
-exceptiongroup==1.1.3
+geopandas==0.13.2
 Fiona==1.9.4.post1
-fsspec==2021.11.1
 GDAL==3.4.1
-geopandas==0.13.2
-greenlet==1.1.2
-HeapDict==1.0.1
-idna==3.4
-importlib-metadata==6.8.0
-importlib-resources==5.4.0
-iniconfig==2.0.0
-Jinja2==3.0.3
-jmespath==0.10.0
-jsonschema==4.2.1
-lark-parser==0.12.0
-locket==0.2.1
-MarkupSafe==2.1.1
-moto==4.2.2
-msgpack==1.0.3
-netCDF4==1.5.8
-numpy==1.21.4
-packaging==21.3
-pandas==1.3.5
-partd==1.2.0
-pluggy==1.3.0
-psutil==5.8.0
-psycopg2==2.9.2
-psycopg2-binary==2.9.2
-py-partiql-parser==0.3.6
-pyarrow==13.0.0
-pycparser==2.21
-pyparsing==3.0.6
-pyrsistent==0.18.0
-pytest==7.4.2
-python-dateutil==2.8.2
-pytz==2021.3
-PyYAML==6.0
-rasterio==1.2.10
-requests==2.31.0
-responses==0.23.3
-s3fs==0.4.2
-s3transfer==0.5.0
-Shapely==1.8.0
-six==1.16.0
-snuggs==1.4.7
-sortedcontainers==2.4.0
-tblib==1.7.0
-tomli==2.0.1
-toolz==0.11.2
-tornado==6.1
-tqdm==4.66.1
-types-PyYAML==6.0.12.11
-urllib3==1.26.7
-werkzeug==2.3.7
-xarray==0.20.2
-xmltodict==0.13.0
-zict==2.0.0
-zipp==3.6.0

dea_conflux/__main__.py

@@ -19,7 +19,6 @@
 import geopandas as gpd
 from datacube.ui import click as ui
 from rasterio.errors import RasterioIOError
-
 import dea_conflux.__version__
 import dea_conflux.db
 import dea_conflux.drill
@@ -176,7 +175,7 @@ def load_and_reproject_shapefile(
     has_s3 = "s3" in gpd.io.file._VALID_URLS
     gpd.io.file._VALID_URLS.discard("s3")
     logger.info(f"Attempting to read {shapefile} to load polgyons.")
-    shapefile = gpd.read_file(shapefile, driver="ESRI Shapefile")
+    shapefile = gpd.read_file(shapefile)
     if has_s3:
         gpd.io.file._VALID_URLS.add("s3")
 
@@ -332,7 +331,6 @@ def run_one(
     Run dea-conflux on one scene.
     """
     logging_setup(verbose)
-
     # Read the plugin as a Python module.
     plugin = run_plugin(plugin)
     logger.info(f"Using plugin {plugin.__file__}")
@@ -703,7 +701,7 @@ def nrt_run_from_queue(
 
             logger.info(f"Processing {id_} ({i + 1}/{len(ids)})")
             logger.info(f"Processing {id_} ({i + 1}/{len(ids)})")
-            
+
             # if id_ is uuid
             if len(id_) != 36:
                 import json
@@ -1384,7 +1382,7 @@ def db_to_csv(
     has_s3 = "s3" in gpd.io.file._VALID_URLS
     gpd.io.file._VALID_URLS.discard("s3")
     logger.info(f"Attempting to read {shapefile} to load polgyons.")
-    shapefile = gpd.read_file(shapefile, driver="ESRI Shapefile")
+    shapefile = gpd.read_file(shapefile)
     if has_s3:
         gpd.io.file._VALID_URLS.add("s3")
 

dea_conflux/drill.py

@@ -19,10 +19,10 @@
 import numpy as np
 import pandas as pd
 import rasterio.features
-import shapely.geometry
 import tqdm
 import xarray as xr
 from datacube.utils.geometry import assign_crs
+from shapely import LineString, Polygon, box
 
 from dea_conflux.types import CRS
 
@@ -106,10 +106,10 @@ def _get_directions(og_geom, int_geom):
 
     Arguments
     ---------
-    og_geom : shapely.geometry.Polygon
+    og_geom : Polygon
         Original polygon.
 
-    int_geom : shapely.geometry.Polygon
+    int_geom : Polygon
         Polygon after intersecting with extent.
 
     Returns
@@ -127,7 +127,7 @@ def _get_directions(og_geom, int_geom):
     for line_ in boundary_intersection_lines:
         coords = list(line_.coords)
         for a, b in zip(coords[:-1], coords[1:]):
-            boundary_intersection_lines_.append(shapely.geometry.LineString((a, b)))
+            boundary_intersection_lines_.append(LineString((a, b)))
     boundary_intersection_lines = boundary_intersection_lines_
 
     boundary_directions = set()
@@ -183,7 +183,7 @@ def _get_directions(og_geom, int_geom):
 
 
 def get_intersections(
-    gdf: gpd.GeoDataFrame, extent: shapely.geometry.Polygon
+    gdf: gpd.GeoDataFrame, extent: Polygon
 ) -> gpd.GeoDataFrame:
     """Find which polygons intersect with an extent and in what direction.
 
@@ -192,7 +192,7 @@ def get_intersections(
     gdf : gpd.GeoDataFrame
         Set of polygons.
 
-    extent : shapely.geometry.Polygon
+    extent : Polygon
         Extent polygon to check intersection against.
 
     Returns
@@ -364,7 +364,7 @@ def filter_shapefile_intersections(
     width = right - left
     height = top - bottom
 
-    testbox = shapely.geometry.Polygon(
+    testbox = Polygon(
         [
             (left - width, bottom - height),
             (left - width, top + height),
@@ -388,12 +388,15 @@ def filter_dataset(dss, shapefile, worker_num=1):
     -------
     filtered_datasets: [str]
     """
-    with multiprocessing.Pool(processes=worker_num) as pool:
-        filtered_datasets = list(
-            tqdm.tqdm(pool.imap(partial(polygon_in_dataset, shapefile=shapefile), dss))
-        )
+    if __name__ == "__main__":
+        with multiprocessing.Pool(processes=worker_num) as pool:
+            filtered_datasets = list(
+                tqdm.tqdm(
+                    pool.imap(partial(polygon_in_dataset, shapefile=shapefile), dss)
+                )
+            )
 
-    return [e for e in filtered_datasets if e]
+        return [e for e in filtered_datasets if e]
 
 
 def polygon_in_dataset(ds, shapefile):
@@ -569,7 +572,7 @@ def drill(
         # of the polygons.
         reference_product = reference_dataset.type.name
         geopolygon = datacube.utils.geometry.Geometry(
-            shapely.geometry.box(*shapefile.total_bounds), crs=shapefile.crs
+            box(*shapefile.total_bounds), crs=shapefile.crs
         )
         time_span = (
             reference_dataset.center_time - time_buffer,

dea_conflux/stack.py

@@ -32,6 +32,9 @@
 import dea_conflux.io
 from dea_conflux.db import Engine
 from dea_conflux.io import CSV_EXTENSIONS, PARQUET_EXTENSIONS
+import dea_tools.bandindices
+import dea_tools.datahandling
+import dea_tools.wetlands
 
 logger = logging.getLogger(__name__)
 
@@ -197,8 +200,10 @@ def remove_timeseries_with_duplicated(df: pd.DataFrame) -> pd.DataFrame:
     # Remove entries within 60s, deals with edge cases where duplicates wrap across midnight UTC
 
     if "date" in df.columns and len(df) > 1:
-        df['TIMEDIFF'] = pd.to_datetime(df['date'].shift(-1)) - pd.to_datetime(df['date'])
-        df = df[~(df['TIMEDIFF'] < timedelta(seconds=60))]
+        df["TIMEDIFF"] = pd.to_datetime(df["date"].shift(-1)) - pd.to_datetime(
+            df["date"]
+        )
+        df = df[~(df["TIMEDIFF"] < timedelta(seconds=60))]
         df = df.drop(columns=["TIMEDIFF"])
 
     # Remember to remove the temp column day in result_df
@@ -223,6 +228,7 @@ def load_pq_file(path):
         date = dea_conflux.io.string_to_date(df.attrs["date"])
         date = stack_format_date(date)
         df.loc[:, "date"] = date
+        df.loc[:, "ard_product"] = str(path).split("/")[-1].split("_")[-3]
     return df
 
 
@@ -279,12 +285,20 @@ def save_df_as_csv(single_polygon_df, feature_id, outpath, remove_duplicated_dat
             / single_polygon_df.loc[norm_veg_index, "overall_veg_num"]
             * single_polygon_df.loc[norm_veg_index, "veg_areas"]
         )
+    # cleaning up single_polygon_df output for stack plot generation
+    single_polygon_df_for_stackplot = single_polygon_df[~(single_polygon_df["pc_missing"] > 0.1)]
+    single_polygon_df_for_stackplot = single_polygon_df_for_stackplot.reset_index()
+    single_polygon_df_for_stackplot["date"] = pd.to_datetime(
+        single_polygon_df_for_stackplot["date"]
+    ).dt.tz_localize(None)
+    dea_tools.wetlands.display_wit_stack_with_df(
+        single_polygon_df_for_stackplot, feature_id, feature_id, x_axis_labels="years"
+    )
 
     # remove the temp column
     single_polygon_df.drop(
         ["overall_veg_num", "veg_areas", "index"], axis=1, inplace=True
     )
-
     if not outpath.startswith("s3://"):
         os.makedirs(Path(filename).parent, exist_ok=True)
     with fsspec.open(filename, "w") as f:
@@ -309,9 +323,9 @@ def stack_wit_tooling_to_single_file(
 
     verbose : bool
     """
+
     polygon_df_list = []
     logger.info("Reading...")
-
     # Note: the stack_wit_tooling_to_single_file() input files are CSV file, which generate by save_df_as_csv()
     # then we assume they already had the norm_pv, norm_npv, norm_bs there.
     with tqdm(total=len(paths)) as bar:
@@ -385,7 +399,8 @@ def stack_wit_tooling(
     """
     wit_df_list = []
     logger.info("Reading...")
-
+    # Note: the stack_wit_tooling_to_single_file() input files are CSV file, which generate by save_df_as_csv()
+    # then we assume they already had the norm_pv, norm_npv, norm_bs there.
     with tqdm(total=len(paths)) as bar:
         with concurrent.futures.ThreadPoolExecutor(max_workers=32) as executor:
             wit_df_list = []
@@ -418,7 +433,6 @@ def stack_wit_tooling(
 
     # delete the temp result to release RAM
     del wit_result
-
     with tqdm(total=len(feature_ids)) as bar:
         with concurrent.futures.ThreadPoolExecutor(
             max_workers=multiprocessing.cpu_count()

examples/wit_ls9.conflux.py

@@ -73,7 +73,7 @@ def transform(inputs: xr.Dataset) -> xr.Dataset:
     # not apply poly_raster cause we will do it before summarise
 
     open_water = wo_ds.water & (1 << 7) > 0
-    
+
     # Thresholding
     # set wet pixels where not masked and above threshold of -350
     wet = tcw.where(mask) > -350

requirements.txt

@@ -1,7 +1,9 @@
-boto3==1.20.24
-datacube==1.8.6
-psycopg2-binary==2.9.2
-pyproj==3.2.1
+boto3==1.35.23
+datacube==1.8.19
+psycopg2-binary==2.9.9
+pyproj==3.7.0
 python-geohash==0.8.5
 s3fs==0.4.2
-SQLAlchemy==1.4.27
+SQLAlchemy==1.4.49
+dea-tools==0.3.5
+seaborn==0.13.2

setup.py

@@ -20,10 +20,10 @@
             "pyarrow",
             "fsspec",
             "s3fs",
-            "moto[s3]",
+            "moto[all]",
             "tqdm",
             "SQLAlchemy",
-            # "python-geohash",
+            "python-geohash",
         ],
         entry_points={
             "console_scripts": ["dea-conflux=dea_conflux.__main__:main"],

tests/data/sum_pv_wet.conflux.py

@@ -1,19 +1,21 @@
 import xarray as xr
 
-product_name = 'sum_pv_wet'
-version = '0.0.1'
+product_name = "sum_pv_wet"
+version = "0.0.1"
 
 input_products = {
-    'ga_ls_wo_3': ['water'],
-    'ga_ls_fc_3': ['PV'],
+    "ga_ls_wo_3": ["water"],
+    "ga_ls_fc_3": ["PV"],
 }
 
 
 def transform(inputs: xr.Dataset) -> xr.Dataset:
-    return xr.Dataset({
-        'water': inputs.water == 128,
-        'pv': inputs.PV > 0.5,
-    })
+    return xr.Dataset(
+        {
+            "water": inputs.water == 128,
+            "pv": inputs.PV > 0.5,
+        }
+    )
 
 
 def summarise(inputs: xr.Dataset) -> xr.Dataset: