Update the DSWx-HLS SAS for cal/val release (#21)

* fill DSWx-HLS products metadata fields: `dem_file`, `landcover_file`, and `worldcover_file`

* fix selection of the resampling algorithm when building overviews

* add metadata to RGB and infrared RGB layers

* add extra margin to the DEM to compute the shadow layer; add metadata to the DEM layer

* fix generation of DSWx-HLS layer overviews

* update DEM_MARGIN_IN_PIXELS from 5 to 10

* update dswx_landcover_mask.py to use files in geographic coordinates as reference for defining the bounding box; add option to define manually the bbox

* fix exclusion of external overviews

* update logging messages

* update dswx_landcover_mask.py bbox argument docstrings

* set default for full_log_formatting to False rather than None

* revert changes to full_log_formatting

* save DSWx-HLS DEM layer as float32 rather than uint8

* save cloud/cloud-shadow masked class with value 9 rather than 2 into binary layer BWTR

* log stderr

* update BWTR color table for cloud/cloud-shadow masked class

* save interpreted layers with 2 water classes rather than 4

* refactor code and update output classes to follow DSWx-HLS product specs;

* force metatadata field SENSING_TIME to have a single data-time entry

* remove temporary files

* update Docker files and build script for cal_val delivery

* fix terrain masking for 4 water classes (before collapsing the 4 water classes into 2)

* address codacy issues

* address codacy issues (2)

* update Zenodo URL with new golden dataset for cal/val delivery

* update ctable for CONF layer class 254 (cloud/cloud-shadow mask); update golden dataset DEM

* update Docker image tag from beta to cal_val (delivery)

* add `--network host` to the docker run command

bin/dswx_landcover_mask.py

@@ -12,7 +12,8 @@
 import argparse
 import logging
 import os
-from osgeo import gdal
+from osgeo import gdal, osr
+import numpy as np
 from proteus.dswx_hls import create_landcover_mask, create_logger
 
 logger = logging.getLogger('dswx_hls_landcover_mask')
@@ -25,14 +26,24 @@ def _get_parser():
         formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
     # Inputs
-    parser.add_argument('input_file',
+    parser.add_argument('-g', '--ref',
+                        dest='reference_file',
                         type=str,
-                        help='Input HLS product')
-
-    parser.add_argument('-c',
-                        '--copernicus-landcover-100m',
+                        help='Reference file for geographic grid'
+                              ' (e.g. HLS product layer)')
+
+    parser.add_argument('-b', '--bbox',
+                        type=float,
+                        nargs=4,
+                        dest='bbox',
+                        metavar=('LAT_MIN', 'LAT_MAX', 'LON_MIN', 'LON_MAX'),
+                        help='Defines the spatial region in '
+                             'the format south north west east.')
+
+    # other parameters
+    parser.add_argument('-c', '--copernicus-landcover-100m',
                         '--landcover', '--land-cover',
-                        dest='landcover_file',
+                        dest='copernicus_landcover_file',
                         required=True,
                         type=str,
                         help='Input Copernicus Land Cover'
@@ -46,8 +57,7 @@ def _get_parser():
                         help='Input ESA WorldCover 10m')
 
     # Outputs
-    parser.add_argument('-o',
-                        '--output-file',
+    parser.add_argument('-o', '--output-file',
                         dest='output_file',
                         required=True,
                         type=str,
@@ -78,24 +88,126 @@ def _get_parser():
     return parser
 
 
+def point2epsg(lon, lat):
+    if lon >= 180.0:
+        lon = lon - 360.0
+    if lat >= 60.0:
+        return 3413
+    elif lat <= -60.0:
+        return 3031
+    elif lat > 0:
+        return 32601 + int(np.round((lon + 177) / 6.0))
+    elif lat < 0:
+        return 32701 + int(np.round((lon + 177) / 6.0))
+    raise ValueError(
+        'Could not determine projection for {0},{1}'.format(lat, lon))
+
 def main():
     parser = _get_parser()
 
     args = parser.parse_args()
 
     create_logger(args.log_file)
 
-    logger.info(f'Input file: {args.input_file}')
-    if not os.path.isfile(args.input_file):
-        logger.error(f'ERROR file not found: {args.input_file}')
+    if args.bbox is None and args.reference_file is None:
+        logger.error('ERROR please, provide either a '
+                     'reference file or a bounding box')
         return
-    layer_gdal_dataset = gdal.Open(args.input_file, gdal.GA_ReadOnly)
-    if layer_gdal_dataset is None:
-        logger.error(f'ERROR invalid input HLS file: {args.input_file}')
-    geotransform = layer_gdal_dataset.GetGeoTransform()
-    projection = layer_gdal_dataset.GetProjection()
-    length = layer_gdal_dataset.RasterYSize
-    width = layer_gdal_dataset.RasterXSize
+
+    if args.bbox is not None:
+        flag_is_geographic = True
+        lat_min, lat_max, lon_min, lon_max = args.bbox
+        # make sure lat_max > lat_min (since dy < 0)
+        if lat_max < lat_min:
+            lat_min, lat_max = lat_max, lat_min
+
+    geographic_srs = None
+    if args.reference_file:
+        print(f'Reference file: {args.reference_file}')
+
+        if not os.path.isfile(args.reference_file):
+            logger.error(f'ERROR file not found: {args.reference_file}')
+            return
+        layer_gdal_dataset = gdal.Open(args.reference_file, gdal.GA_ReadOnly)
+        if layer_gdal_dataset is None:
+            logger.error(f'ERROR invalid file: {args.reference_file}')
+
+        # read reference image geolocation
+        geotransform = layer_gdal_dataset.GetGeoTransform()
+        projection = layer_gdal_dataset.GetProjection()
+        length = layer_gdal_dataset.RasterYSize
+        width = layer_gdal_dataset.RasterXSize
+
+        # check if geolocation is in projected coordinates, i.e.,
+        # not geographic
+        projection_srs = osr.SpatialReference(wkt=projection)
+
+        flag_is_geographic = projection_srs.IsGeographic()
+        if flag_is_geographic:
+            geographic_srs = projection_srs
+            print('Reference file is provided in geographic'
+                  ' coordinates.')
+            lat_max = geotransform[3]
+            lat_min = lat_max + geotransform[5] * length
+            lon_min = geotransform[0]
+            lon_max = lon_min + geotransform[1] * width
+        else:
+            print('Reference file is NOT provided in geographic'
+                        ' coordinates')
+
+    if flag_is_geographic:
+        mean_y = (lat_max + lat_min) / 2.0
+        mean_x = (lon_min + lon_max) / 2.0
+        epsg = point2epsg(mean_x, mean_y)
+
+        print('Converting geographic coordinates to UTM:')
+        print(f'    lat_min: {lat_min}, lat_max: {lat_max}')
+        print(f'    lon_min: {lon_min}, lon_max: {lon_max}')
+        print(f'    EPSG code: {epsg}')
+
+        if geographic_srs is None:
+            geographic_srs = osr.SpatialReference()
+            geographic_srs.SetWellKnownGeogCS("WGS84")
+
+        utm_srs = osr.SpatialReference()
+        utm_srs.ImportFromEPSG(epsg)
+
+        # create transformation of coordinates from geographic (lat/lon)
+        # to UTM
+        transformation = osr.CoordinateTransformation(geographic_srs, utm_srs)
+        y_min = None
+        y_max = None
+        x_min = None
+        x_max = None
+        for lat in [lat_max, lat_min]:
+            for lon in [lon_min, lon_max]:
+                x, y, _ = transformation.TransformPoint(lat, lon, 0)
+                if y_min is None or y_min > y:
+                    y_min = y
+                if y_max is None or y_max < y:
+                    y_max = y
+                if x_min is None or x_min > x:
+                    x_min = x
+                if x_max is None or x_max < x:
+                    x_max = x
+
+        print(f'    y_min: {y_min}, y_max: {y_max}')
+        print(f'    x_min: {x_min}, x_max: {x_max}')
+
+        # land cover map step is 30m meters
+        dx = 30
+        dy = -30
+        geotransform = [x_min, dx, 0, y_max, 0, dy]
+
+        width = int(np.ceil((x_max - x_min) / dx))
+        length = int(np.ceil((y_min - y_max) / dy))
+        projection = utm_srs.ExportToProj4()
+        projection = projection.strip()
+
+        print(f'    width: {width}')
+        print(f'    length: {length}')
+        print(f'    geotransform: {geotransform}')
+        print(f'    projection: {projection}')
 
     create_landcover_mask(args.copernicus_landcover_file,
                           args.worldcover_file, args.output_file,

build_docker_image.sh

@@ -1,7 +1,7 @@
 #!/bin/bash
 
 IMAGE=opera/proteus
-t=beta
+t=cal_val
 echo "IMAGE is $IMAGE:$t"
 
 # fail on any non-zero exit codes
@@ -13,7 +13,7 @@ python3 setup.py sdist
 docker build --rm --force-rm --network=host -t ${IMAGE}:$t -f docker/Dockerfile .
 
  # run tests
-docker run --rm -u "$(id -u):$(id -g)" -v "$PWD:/mnt" -w /mnt -it "${IMAGE}:$t" pytest /mnt/tests/
+docker run --rm -u "$(id -u):$(id -g)" -v "$PWD:/mnt" -w /mnt -it --network host "${IMAGE}:$t" pytest /mnt/tests/
 
 # create image tar
 docker save opera/proteus > docker/dockerimg_proteus_$t.tar

docker/Dockerfile

@@ -25,7 +25,7 @@ RUN set -ex \
  && conda update --all --yes \
  && conda install --yes --file /tmp/requirements.txt \
  && conda install --yes --channel conda-forge --file /tmp/requirements_forge.txt \
- && conda clean -tipsy \
+ && conda clean -tip \
  && rm -rf /opt/conda/pkgs \
  && rm /tmp/requirements.txt \
  && rm /tmp/requirements_forge.txt

load_docker_tar.sh

@@ -1,3 +1,3 @@
 #!/bin/bash
 
-docker load -i docker/dockerimg_proteus_beta.tar
+docker load -i docker/dockerimg_proteus_cal_val.tar

src/proteus/core.py

@@ -40,34 +40,33 @@ def save_as_cog(filename, scratch_dir = '.', logger = None,
     else:
         resamp_algorithm = 'CUBICSPLINE'
 
-    gdal_ds.BuildOverviews('CUBICSPLINE', overviews_list,
+    gdal_ds.BuildOverviews(resamp_algorithm, overviews_list,
                            gdal.TermProgress_nocb)
 
     del gdal_ds  # close the dataset (Python object and pointers)
     external_overview_file = filename + '.ovr'
     if os.path.isfile(external_overview_file):
-        os.path.remove(external_overview_file)
+        os.remove(external_overview_file)
 
     logger.info('COG step 2: save as COG')
     temp_file = tempfile.NamedTemporaryFile(
                     dir=scratch_dir, suffix='.tif').name
 
     tile_size = 256
-    ovr_tile_size = tile_size
-    gdal_translate_options = [
-        'TILED=YES',
-        f'BLOCKXSIZE={tile_size}',
-        f'BLOCKYSIZE={tile_size}',
-        f'GDAL_TIFF_OVR_BLOCKSIZE={ovr_tile_size}'
-        'COPY_SRC_OVERVIEWS=YES'] 
+    # ovr_tile_size = tile_size
+    gdal_translate_options = ['TILED=YES',
+                              f'BLOCKXSIZE={tile_size}',
+                              f'BLOCKYSIZE={tile_size}',
+                              # f'GDAL_TIFF_OVR_BLOCKSIZE={ovr_tile_size}'
+                              'COPY_SRC_OVERVIEWS=YES'] 
 
     if flag_compress:
         gdal_translate_options += ['COMPRESS=DEFLATE']
 
     if is_integer:
         gdal_translate_options += ['PREDICTOR=2']
     else:
-        gdal_translate_options = ['PREDICTOR=3']
+        gdal_translate_options += ['PREDICTOR=3']
 
     gdal.Translate(temp_file, filename,
                    creationOptions=gdal_translate_options)
@@ -102,7 +101,6 @@ def get_geographic_boundaries_from_mgrs_tile(mgrs_tile_name, verbose=False):
 
     # create UTM spatial reference
     utm_coordinate_system = osr.SpatialReference()
-    utm_coordinate_system.SetWellKnownGeogCS("WGS84")
     utm_coordinate_system.SetUTM(utm_zone, is_northern)
 
     # create geographic (lat/lon) spatial reference