Improvements to coastal masking, add mask modifications

coastlines/continental.py

@@ -136,6 +136,8 @@ def continental_cli(
  baseline_year,
  include_styles,
 ):
+
+ print(hotspots)
 
  #################
  # Merge vectors #

coastlines/vector.py

@@ -483,7 +483,7 @@ def contours_preprocess(
  gapfill_ds,
  water_index,
  index_threshold,
- buffer_pixels=33,
+ buffer_pixels=50,
  mask_temporal=True,
  mask_modifications=None,
  debug=False,
@@ -577,27 +577,34 @@ def contours_preprocess(
  temporal_mask = temporal_masking(thresholded_ds == 1)
  thresholded_ds = thresholded_ds.where(temporal_mask)
 
- # Create all time layer by identifying pixels that are land in at
- # least 15% of valid observations; this is used to generate an
- # all-of-time buffered coastal study area that constrains the Coastlines
- # analysis to the coastal zone
- all_time = thresholded_ds.mean(dim="year") > 0.5 # 0.15
-
- # Remove narrow river and stream features from all time layer using
- # the `black_tophat` transform. So narrow channels between islands
- # and the mainland aren't mistaken for rivers, first apply `sieve`
- # to any connected land pixels (i.e. islands) smaller than 5 km^2
+ # Create all time layers by identifying pixels that are land in at
+ # least 15% and 50% of valid observations; the 15% layer is used to
+ # identify pixels that contain land for even a small period of time;
+ # this is used for producing an all-time buffered "coastal" study area.
+ # The 50% layer is used to more conservatively identify pixels that
+ # are land for a majority of years; this is used for extracting
+ # rivers.
+ all_time_15 = thresholded_ds.mean(dim="year") > 0.15
+ all_time_50 = thresholded_ds.mean(dim="year") > 0.5
+
+ # Identify narrow river and stream features using the `black_tophat`
+ # transform. To avoid narrow channels between islands and the
+ # mainland being mistaken for rivers, first apply `sieve` to any
+ # connected land pixels (i.e. islands) smaller than 5 km^2.
+ # Use a disk of size 8 to identify any rivers/streams smaller than
+ # approximately 240 m (e.g. 8 * 30 m = 240 m).
  island_size = int(5000000 / (30 * 30)) # 5 km^2
- sieved = xr.apply_ufunc(sieve, all_time.astype(np.int16), island_size)
- rivers = xr.apply_ufunc(black_tophat, (sieved & all_time), disk(5)) == 1
- all_time_clean = all_time.where(~rivers, True)
+ sieved = xr.apply_ufunc(sieve, all_time_50.astype(np.int16), island_size)
+ rivers = xr.apply_ufunc(black_tophat, (sieved & all_time_50), disk(8)) == 1
 
  # Create a river mask by eroding the all time layer to clip out river
- # and estuary mouths, then expanding river features to include stream
- # banks and account for migrating rivers
- all_time_eroded = xr.apply_ufunc(binary_erosion, all_time_clean, disk(10))
- rivers = rivers.where(all_time_eroded, False)
- river_mask = ~xr.apply_ufunc(binary_dilation, rivers, disk(5))
+ # mouths, then expanding river features to include stream banks and
+ # account for migrating rivers
+ river_mouth_mask = xr.apply_ufunc(
+ binary_erosion, all_time_50.where(~rivers, True), disk(12)
+ )
+ rivers = rivers.where(river_mouth_mask, False)
+ river_mask = ~xr.apply_ufunc(binary_dilation, rivers, disk(4))
 
  # Load Geodata 100K coastal layer to use to separate ocean waters from
  # other inland waters. This product has values of 0 for ocean waters,
@@ -613,24 +620,46 @@ def contours_preprocess(
  ocean_da = xr.apply_ufunc(binary_erosion, geodata_da == 0, disk(10))
 
  # Use all time and Geodata 100K data to produce the buffered coastal
- # study area. Morphological closing helps to "close" the entrances of
- # estuaries and rivers, removing them from the analysis.
+ # study area. The output has values of 0 representing non-coastal
+ # "ocean", values of 1 representing "coastal", and values of 2
+ # representing non-coastal "inland" pixels.
  coastal_mask = coastal_masking(
- ds=all_time_clean, ocean_da=ocean_da, buffer=buffer_pixels
+ ds=all_time_15.where(~rivers, True), ocean_da=ocean_da, buffer=buffer_pixels
  )
 
- # The output of `coastal_masking` contains values of 2 that represent
- # pixels inland of the coastal buffer. Use this to create
- # a mask of inland regions:
- inland_mask = coastal_mask != 2
+ # Add rivers as "inland" pixels in the coastal mask
+ coastal_mask = xr.where(river_mask, coastal_mask, 2)
+
+ # Optionally modify the coastal mask using manually supplied
+ # polygons to add missing areas of shoreline, or remove unwanted
+ # areas from the mask.
+ if mask_modifications is not None:
+
+ # Only proceed if there are polygons available
+ if len(mask_modifications.index) > 0:
+
+ # Convert type column to integer, with 1 representing pixels
+ # to add to the coastal mask (by setting them as "coastal"
+ # pixels, and 2 representing pixels to remove from the mask
+ # (by setting them as "inland data")
+ mask_modifications = mask_modifications.replace({"add": 1, "remove": 2})
+
+ # Rasterise polygons into extent of satellite data
+ modifications_da = xr_rasterize(
+ mask_modifications, da=yearly_ds, attribute_col="type"
+ )
+
+ # Where `modifications_da` has a value other than 0,
+ # replace values from `coastal_mask` with `modifications_da`
+ coastal_mask = coastal_mask.where(modifications_da == 0, modifications_da)
 
  # Generate individual annual masks by selecting only water pixels that
  # are directly connected to the ocean in each yearly timestep
  annual_mask = (
  # Treat both 1s and NaN pixels as land (i.e. True)
  (thresholded_ds != 0)
- # Mask out rivers and inland regions
- .where(river_mask & inland_mask)
+ # Mask out inland regions (i.e. values of 2 in `coastal_mask`)
+ .where(coastal_mask != 2)
  # Keep pixels directly connected to ocean in each timestep
  .groupby("year").map(
  func=ocean_masking,
@@ -641,9 +670,9 @@ def contours_preprocess(
  )
 
  # Finally, apply temporal and annual masks to our surface water
- # index data, then clip to the all-of-time coastal study area
+ # index data, then clip to "coastal" pixels in `coastal_mask`
  masked_ds = combined_ds[water_index].where(
- temporal_mask & annual_mask & coastal_mask
+ temporal_mask & annual_mask & (coastal_mask == 1)
  )
 
  # Generate annual vector polygon masks containing information
@@ -655,14 +684,14 @@ def contours_preprocess(
  return (
  masked_ds,
  certainty_masks,
- all_time,
- all_time_clean,
+ all_time_15,
+ all_time_50,
  river_mask,
  ocean_da,
- coastal_mask,
- inland_mask,
  thresholded_ds,
+ temporal_mask,
  annual_mask,
+ coastal_mask,
  )
 
  else:
@@ -1391,10 +1420,10 @@ def generate_vectors(
  gridcell_gdf.index = gridcell_gdf.index.astype(int).astype(str)
  gridcell_gdf = gridcell_gdf.loc[[str(study_area)]]
 
- #  # Coastal mask modifications
- # modifications_gdf = gpd.read_file(
- #  config["Input files"]["modifications_path"], bbox=bbox
- # ).to_crs(str(yearly_ds.odc.crs))
+ # Coastal mask modifications
+ modifications_gdf = gpd.read_file(
+ config["Input files"]["modifications_path"], bbox=bbox
+ ).to_crs(str(yearly_ds.odc.crs))
 
  # Geomorphology dataset
  geomorphology_gdf = gpd.read_file(
@@ -1417,7 +1446,7 @@ def generate_vectors(
  water_index,
  index_threshold,
  buffer_pixels=33,
- mask_modifications=None,
+ mask_modifications=modifications_gdf,
  )
 
  # Extract annual shorelines

configs/dea_coastlines_config.yaml

@@ -3,6 +3,6 @@ Virtual product:
  virtual_product_name: ls_nbart_ndwi
 Input files:
  grid_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/albers_grids/ga_summary_grid_c3_48km_coastal_clipped.geojson
- modifications_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/estuary_mask_modifications.geojson
+ modifications_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/modifications_deacoastlines.geojson
  geomorphology_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/Smartline.gpkg
  region_attributes_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/Primary_compartments.geojson

configs/dea_coastlines_config_development.yaml

@@ -3,6 +3,6 @@ Virtual product:
  virtual_product_name: ls_nbart_ndwi
 Input files:
  grid_path: data/raw/coastlines_development.geojson
- modifications_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/estuary_mask_modifications.geojson
+ modifications_path: data/raw/modifications_deacoastlines.geojson
  geomorphology_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/Smartline.gpkg
  region_attributes_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/Primary_compartments.geojson

configs/dea_coastlines_config_tests.yaml

@@ -3,6 +3,6 @@ Virtual product:
  virtual_product_name: ls_nbart_ndwi
 Input files:
  grid_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/albers_grids/ga_summary_grid_tests.geojson
- modifications_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/estuary_mask_modifications.geojson
+ modifications_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/modifications_deacoastlines.geojson
  geomorphology_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/Smartline.gpkg
  region_attributes_path: https://dea-public-data.s3.ap-southeast-2.amazonaws.com/derivative/dea_coastlines/supplementary/Primary_compartments.geojson