Tests and fleshing it all out

Author: omad

Dockerfile

@@ -1,5 +1,5 @@
 
-# Create image with docker buildx build -t datacube-vectoriser .
+# Create image with docker buildx build -t dea-vectoriser .
 
 # Thankyou https://uwekorn.com/2021/03/01/deploying-conda-environments-in-docker-how-to-do-it-right.html
 
@@ -21,4 +21,4 @@ RUN --mount=type=cache,target=/opt/conda/pkgs micromamba install -y -n base -f /
     find /opt/conda/ -name '__pycache__' -type d -exec rm -rf '{}' '+'
 
 
-COPY . .
+COPY . .

datacube_vectoriser/utils.py

@@ -1,14 +1,17 @@
-import json
-import os
 from concurrent import futures
 
 import boto3
-from toolz import dicttoolz
+import json
+import os
+from pathlib import PurePosixPath
+from toolz import dicttoolz, get_in
+from typing import Tuple, Optional
+from urllib.parse import urlparse
 
 
 def _stac_to_sns(sns_arn, stac):
     """
-    Publish our STAC document to an SNS
+    Publish a STAC document to an SNS
     """
     bbox = stac["bbox"]
 
@@ -41,6 +44,7 @@ def _stac_to_sns(sns_arn, stac):
 
 
 def upload_directory(directory, bucket, prefix, boto3_session: boto3.Session = None):
+    """Recursively upload a directory to an s3 bucket"""
     if boto3_session is None:
         boto3_session = boto3.Session()
     s3 = boto3_session.client("s3")
@@ -71,7 +75,9 @@ def upload_file(filename):
             except Exception as e:
                 print("Exception {} encountered while uploading file {}".format(e, upload_task[task]))
 
+
 def receive_messages(queue_url):
+    """Yield SQS Messages until the queue is empty"""
     sqs = boto3.resource('sqs')
     queue = sqs.Queue(queue_url)
 
@@ -80,13 +86,32 @@ def receive_messages(queue_url):
 
     while len(messages) > 0:
         for message in messages:
-            # body = json.loads(message.body)
             yield message
 
-            # message.delete()
-
         messages = queue.receive_messages(MaxNumberOfMessages=1, )
 
+
+def geotiff_url_from_stac(stac_document) -> Optional[str]:
+    return get_in(['assets', 'water', 'href'], stac_document)
+
+
+def output_name_from_url(src_url, file_suffix) -> Tuple[PurePosixPath, str]:
+    """Derive the output directory structure and filename from the input URL"""
+    o = urlparse(src_url)
+    path = PurePosixPath(o.path)
+
+    relative_path = PurePosixPath(*path.parts[-6:-1])
+
+    filename = path.with_suffix(file_suffix).name
+
+    # parts
+    # Out[6]: ['097', '075', '1998', '08', '17']
+    # filename
+    # Out[7]: 'ga_ls_wo_3_097075_1998-08-17_final_water.tif'
+
+    return relative_path, filename
+
+
 def chain_funcs(arg, *funcs):
     result = arg
     for f in funcs:

datacube_vectoriser/vector_wos.py

@@ -6,81 +6,83 @@
 #     - C) conduct 1 pixle buffer of no-data class? (unsure if should be latter in workflow)
 # 3. vectorise
 # 4. simplify shapes to remove complixity
-# 5. join both data types back together as one Geopandas Geodataframe (container for sapely objects with projection imformation)
+# 5. join both data types back together as one Geopandas Geodataframe (container for shapely objects with projection
+#    information)
 # 6. export an a single shapefile with attributes intact.
 
-import fiona
-# Derived from https://github.com/GeoscienceAustralia/dea-notebooks/blob/KooieCate/vector_WOs_draft4.py
 import geopandas
 import geopandas as gp
 import pandas as pd
 import rasterio.features
 import xarray as xr
 from fiona.crs import from_epsg
 from scipy import ndimage
-from shapely.geometry import shape, mapping
+from shapely.geometry import shape
+# Derived from https://github.com/GeoscienceAustralia/dea-notebooks/blob/KooieCate/vector_WOs_draft4.py
+from typing import Tuple
 
 
-def load_data(url):
-    # Open geotiff and reformat to Xarray DataArray
-    geotiff_wos = xr.open_rasterio(url)  # 'ga_s2am_wo_0-0-1_49JGN_2021-02-08_nrt_water.tif')
+def load_wos_data(url) -> xr.Dataset:
+    """Open a GeoTIFF info an in memory DataArray """
+    geotiff_wos = xr.open_rasterio(url)
     wos_dataset = geotiff_wos.to_dataset('band')
     wos_dataset = wos_dataset.rename({1: 'wo'})
     return wos_dataset
 
 
-def generate_raster_layers(wos_dataset):
+def generate_raster_layers(wos_dataset: xr.Dataset) -> Tuple[xr.DataArray, xr.DataArray]:
     # Defining the three 'classes':
     # a) Water: where water is observed. Bit value 128
-    # b) unspoken 'dry'. this is not vectorised and is left and transparent layer. bit values: 1 (no data) 2 (Contiguity)
-    # c) Not_analysed: every masking applied to the data except terrain shadow. bit values: composed of Everyting else,
+    # b) unspoken 'dry'. this is not vectorised and is left as a transparent layer. bit values: 1 (no data) 2 (
+    # Contiguity)
+    # c) Not_analysed: every masking applied to the data except terrain shadow. bit values: composed of everything else,
     # 1 create binary arrays for two classes of interest
     water_vals = (wos_dataset.wo == 128)  # water only has 128 water observations
-    # here we used reversed logic to turn all pixles that should be 'not analysed' to a value of 3. is is easier to list the 4 classes that are passed to the unlabled 'dry' class
+    # here we used reversed logic to turn all pixels that should be 'not analysed' to a value of 3. is is easier to
+    # list the 4 classes that are passed to the unlabled 'dry' class
     not_analysed = wos_dataset.wo.where(((wos_dataset.wo == 0) | (wos_dataset.wo == 1) | (wos_dataset.wo == 8)
                                          | (wos_dataset.wo == 2) | (wos_dataset.wo == 128) | (wos_dataset.wo == 130) | (
                                                  wos_dataset.wo == 142)), 3)
-    not_analysed = not_analysed.where((not_analysed == 3), 0)  # now keep the 3 values and make everyting else 0
-    # 2 conduct binary errosion and closing to remove single pixles
+    not_analysed = not_analysed.where((not_analysed == 3), 0)  # now keep the 3 values and make everything else 0
+    # 2 conduct binary erosion and closing to remove single pixels
     erroded_water = xr.DataArray(ndimage.binary_erosion(water_vals, iterations=2).astype(water_vals.dtype),
                                  coords=water_vals.coords)
     erroded_not_analysed = xr.DataArray(ndimage.binary_erosion(not_analysed, iterations=2).astype(not_analysed.dtype),
                                         coords=not_analysed.coords)
-    # dialating cloud 3 times after erroding 2, to create small overlap and iliminate gaps in data
+    # dilating cloud 3 times after eroding 2, to create small overlap and illuminate gaps in data
     dilated_water = xr.DataArray(ndimage.binary_dilation(erroded_water, iterations=3).astype(water_vals.dtype),
                                  coords=water_vals.coords)
     dilated_not_analysed = xr.DataArray(
-        ndimage.binary_dilation(erroded_not_analysed, iterations=(3)).astype(not_analysed.dtype),
+        ndimage.binary_dilation(erroded_not_analysed, iterations=3).astype(not_analysed.dtype),
         coords=not_analysed.coords)
 
     return dilated_water, dilated_not_analysed
 
 
-def vectorise_data(xarrayDataArray, transform, crs, label='Label'):
+def vectorise_data(data_array: xr.DataArray, transform, crs, label='Label'):
     """this module takes an Xarray DataArray and vectorises it as shapely geometries in a Geopandas Geodataframe
 
     Input
-    xarrayDataArray: an Xarray DataArray with boolean values (1,0) with 1 or True equal to the areas that will be turned into vectors
-    Label: default 'Label', String, the data label that will be added to each geometry in geodataframe
+    data_array: a DataArray with boolean values (1,0) with 1 or True equal to the areas that will be turned
+                into vectors
+    label: default 'Label', String, the data label that will be added to each geometry in geodataframe
 
     output
-    Geodataframe containing shapely geometies with data type lable in a series called attribute"""
+    Geodataframe containing shapely geometries with data type label in a series called attribute"""
 
     vector = rasterio.features.shapes(
-        xarrayDataArray.data.astype('float32'),
-        mask=xarrayDataArray.data.astype('float32') == 1,  # this defines which part of array becomes polygons
+        data_array.data.astype('float32'),
+        mask=data_array.data.astype('float32') == 1,  # this defines which part of array becomes polygons
         transform=transform)
 
-    # rasterio.features.shapes outputs tupples. we only want the polygon coordinate portions of the tupples
-    vectored_data = list(vector)  # put tupple output in list
+    # rasterio.features.shapes outputs tuples. we only want the polygon coordinate portions of the tuples
+    vectored_data = list(vector)  # put tuple output in list
 
     # Extract the polygon coordinates from the list
     polygons = [polygon for polygon, value in vectored_data]
-    # create empty list for lables
-    labels = []
-    # put in labels
-    for i in polygons:
-        labels.append(label)  # create a list with the data label type
+
+    # create a list with the data label type
+    labels = [label for _ in polygons]
 
     # Convert polygon coordinates into polygon shapes
     polygons = [shape(polygon) for polygon in polygons]
@@ -93,7 +95,7 @@ def vectorise_data(xarrayDataArray, transform, crs, label='Label'):
 
 
 def vectorise_wos_from_url(url) -> geopandas.GeoDataFrame:
-    raster = load_data(url)
+    raster = load_wos_data(url)
 
     dataset_crs = from_epsg(raster.crs[11:])
     dataset_transform = raster.transform
@@ -102,26 +104,25 @@ def vectorise_wos_from_url(url) -> geopandas.GeoDataFrame:
     dilated_water, dilated_not_analysed = generate_raster_layers(raster)
 
     # vectorise the arrays
-
     notAnalysedGPD = vectorise_data(dilated_not_analysed, dataset_transform, dataset_crs, label='Not_analysed')
 
     WaterGPD = vectorise_data(dilated_water, dataset_transform, dataset_crs, label='Water')
 
     # Simplify
 
-    # Run simplification with 15 tollerance
-    simplifyed_water = WaterGPD.simplify(10)
+    # Run simplification with 15 tolerance
+    simplified_water = WaterGPD.simplify(10)
 
-    simplifyed_notAnalysed = notAnalysedGPD.simplify(15)
+    simplified_not_analysed = notAnalysedGPD.simplify(15)
 
     # Put simplified shapes in a dataframe
-    simple_waterGPD = gp.GeoDataFrame(geometry=simplifyed_water,
+    simple_waterGPD = gp.GeoDataFrame(geometry=simplified_water,
                                       crs=dataset_crs)
 
-    simple_notAnalysedGPD = gp.GeoDataFrame(geometry=simplifyed_notAnalysed,
+    simple_notAnalysedGPD = gp.GeoDataFrame(geometry=simplified_not_analysed,
                                             crs=dataset_crs)
 
-    # add attribute lables back in
+    # add attribute labels back in
     simple_waterGPD['attribute'] = WaterGPD['attribute']
 
     simple_notAnalysedGPD['attribute'] = notAnalysedGPD['attribute']
@@ -132,38 +133,7 @@ def vectorise_wos_from_url(url) -> geopandas.GeoDataFrame:
 
     # 6 Join together and save to file
 
-    All_classes = gp.GeoDataFrame(pd.concat([simple_waterGPD, simple_notAnalysedGPD], ignore_index=True),
+    all_classes = gp.GeoDataFrame(pd.concat([simple_waterGPD, simple_notAnalysedGPD], ignore_index=True),
                                   crs=simple_notAnalysedGPD.crs)
 
-    return All_classes
-
-
-def save_to_file(gpd_dataframe):
-    # +
-    # #define output file name to save vectors as
-    outFile = 'Test_simplify/WO_vectors_test_clean'
-
-    # # Save the polygons to a shapefile
-    schema = {
-        'geometry': 'Polygon',
-        'properties': {
-            'attribute': 'str'
-        }
-    }
-
-    # # Generate our dynamic filename
-    FileName = f'{outFile}.shp'
-
-    # #create file and save
-    with fiona.open(FileName,
-                    "w",
-                    crs=from_epsg(3577),
-                    driver='ESRI Shapefile',
-                    schema=schema) as output:
-        for ix, poly in gpd_dataframe.iterrows():
-            output.write(({
-                'properties': {
-                    'attribute': poly['attribute']
-                },
-                'geometry': mapping(shape(poly['geometry']))
-            }))
+    return all_classes

environment.yaml

@@ -3,7 +3,6 @@ channels:
   - conda-forge
 dependencies:
   - python=3.9
-  - datacube
   - boto3
   - scipy
   - rasterio
@@ -15,4 +14,4 @@ dependencies:
   - xarray
   - toolz
   - pip:
-      - boto3_fixtures
+      - boto3_fixtures

setup.cfg

@@ -0,0 +1,2 @@
+[flake8]
+max-line-length = 120

tests/conftest.py

@@ -1,11 +1,32 @@
-import fixtures
 import boto3_fixtures as b3f
+import pytest
+from pathlib import Path
+from urllib.request import urlretrieve
+
+import fixtures
 
 aws = b3f.contrib.pytest.moto_fixture(
-  services=["sqs", "s3", "sns"],
-  scope="class",
+    services=["sqs", "s3", "sns"],
+    scope="class",
 )
 
 sqs = b3f.contrib.pytest.service_fixture("sqs", scope="class", queues=fixtures.SQS)
 s3 = b3f.contrib.pytest.service_fixture("s3", scope="class", buckets=fixtures.S3)
 sns = b3f.contrib.pytest.service_fixture("sns", scope="class", topics=fixtures.SNS)
+
+SAMPLE_DATA = [
+    ("sample_raster.tif", "https://data.dea.ga.gov.au/derivative/ga_ls_wo_3/1-6-0/097/075/1998/08/17/"
+                          "ga_ls_wo_3_097075_1998-08-17_final_water.tif"),
+    ("sample_stac.json", "https://data.dea.ga.gov.au/derivative/ga_ls_wo_3/1-6-0/097/075/1998/08/17/"
+                         "ga_ls_wo_3_097075_1998-08-17_final.stac-item.json")
+]
+
+
+@pytest.fixture
+def sample_data(pytestconfig):
+    data_dir = Path(pytestconfig.cache.makedir('vect_data'))
+    for dest, src in SAMPLE_DATA:
+        destination_file = data_dir / dest
+        if not destination_file.exists():
+            urlretrieve(src, destination_file)
+    return data_dir

tests/fixtures.py

@@ -13,4 +13,4 @@
             "DisplayName": "YourSystemIsOnFireTopic",
         },
     }
-]
+]