Merge branch 'master' into fix/rde/sentinel-2-z-values

README.md

@@ -98,11 +98,11 @@ Once the container is built, you can run the `scripts/` scripts inside a docker
 
 ### Using virtualenv
 
-If not using docker, itt's recommended to use [virtualenv](https://virtualenv.pypa.io/en/latest/index.html) to keep isolate the python environment used to develop stactools. See virtualenv documentation for more detailed information, but as a shortcut here's some quick steps:
+If not using docker, it's recommended to use [virtualenv](https://virtualenv.pypa.io/en/latest/index.html) to keep isolate the python environment used to develop stactools. See virtualenv documentation for more detailed information, but as a shortcut here's some quick steps:
 
 - Make sure [virtualenv](https://virtualenv.pypa.io/en/latest/installation.html) is installed
 - Run `virtualenv venv`
-- Activate the virtualenv with `source venv/bin/active`
+- Activate the virtualenv with `source venv/bin/activate`
 
 #### Installing development requirements
 

requirements-dev.txt

@@ -10,3 +10,4 @@ yapf==0.30.*
 ipython==7.16.1
 nbsphinx==0.7.1
 coverage==5.2.*
+

stactools_aster/stactools/aster/cog.py

@@ -10,6 +10,7 @@
 
 from stactools.core.projection import reproject_geom
 from stactools.core.utils.subprocess import call
+from stactools.core.utils.convert import cogify
 from stactools.aster.utils import AsterSceneId
 from stactools.aster.xml_metadata import XmlMetadata
 
@@ -41,13 +42,6 @@ def set_band_names(href: str, band_names: List[str]) -> None:
         ds.descriptions = band_names
 
 
-def cogify(input_path, output_path):
-    call([
-        'gdal_translate', '-of', 'COG', '-co', 'predictor=2', '-co',
-        'compress=deflate', input_path, output_path
-    ])
-
-
 def _create_cog_for_sensor(sensor: str, file_prefix: str, tmp_dir: str,
                            output_dir: str, bounds: List[float], crs: str,
                            subdataset_info: List[Tuple[Any, int]]) -> str:
@@ -70,7 +64,7 @@ def _create_cog_for_sensor(sensor: str, file_prefix: str, tmp_dir: str,
 
     set_band_names(merged_path, band_names)
 
-    cogify(merged_path, sensor_cog_href)
+    cogify(merged_path, sensor_cog_href, extra_args=["-co", "predictor=2"])
 
     return sensor_cog_href
 

stactools_cgls_lc100/stactools/cgls_lc100/cog.py

@@ -1,30 +1,15 @@
 import logging
 import os
-from subprocess import Popen, PIPE, STDOUT
 
 import pystac
 
+from stactools.core.utils.convert import cogify
 from stactools.cgls_lc100.constants import (ITEM_COG_IMAGE_NAME,
                                             ITEM_TIF_IMAGE_NAME)
 
 logger = logging.getLogger(__name__)
 
 
-def call(command):
-    def log_subprocess_output(pipe):
-        for line in iter(pipe.readline, b''):  # b'\n'-separated lines
-            logger.info(line.decode("utf-8").strip('\n'))
-
-    process = Popen(command, stdout=PIPE, stderr=STDOUT)
-    with process.stdout:
-        log_subprocess_output(process.stdout)
-    return process.wait()  # 0 means success
-
-
-def cogify(input_path, output_path):
-    call(['gdal_translate', '-of', 'COG', '-co', input_path, output_path])
-
-
 def _create_cog(item, cog_directory):
     tif_asset = item.assets.get(ITEM_TIF_IMAGE_NAME)
     cogify(tif_asset.href, cog_directory)

stactools_core/stactools/core/utils/convert.py

@@ -0,0 +1,20 @@
+from typing import List
+
+from stactools.core.utils.subprocess import call
+
+DEFAULT_COGIFY_ARGS = ["-co", "compress=deflate"]
+
+
+def cogify(infile: str,
+           outfile: str,
+           args: List[str] = None,
+           extra_args: List[str] = None):
+    """Creates a COG from a GDAL-readable file."""
+    if args is None:
+        args = DEFAULT_COGIFY_ARGS[:]
+    args = ["gdal_translate", "-of", "COG"] + args
+    if extra_args:
+        args.extend(extra_args)
+    args.append(infile)
+    args.append(outfile)
+    return call(args)

stactools_corine/stactools/corine/cog.py

@@ -1,33 +1,15 @@
 import logging
 import os
-from subprocess import Popen, PIPE, STDOUT
 
 import pystac
 
+from stactools.core.utils.convert import cogify
 from stactools.corine.constants import (ITEM_COG_IMAGE_NAME,
                                         ITEM_TIF_IMAGE_NAME)
 
 logger = logging.getLogger(__name__)
 
 
-def call(command):
-    def log_subprocess_output(pipe):
-        for line in iter(pipe.readline, b''):  # b'\n'-separated lines
-            logger.info(line.decode("utf-8").strip('\n'))
-
-    process = Popen(command, stdout=PIPE, stderr=STDOUT)
-    with process.stdout:
-        log_subprocess_output(process.stdout)
-    return process.wait()  # 0 means success
-
-
-def cogify(input_path, output_path):
-    call([
-        'gdal_translate', '-of', 'COG', '-co', 'compress=deflate', input_path,
-        output_path
-    ])
-
-
 def _create_cog(item, cog_directory):
 
     tif_asset = item.assets.get(ITEM_TIF_IMAGE_NAME)

stactools_landsat/stactools/landsat/mtl_metadata.py

@@ -2,6 +2,7 @@
 from typing import Any, Dict, List, Optional
 
 from pystac.utils import str_to_datetime
+from pyproj import Geod
 
 from stactools.core.utils import map_opt
 from stactools.core.projection import transform_from_bbox
@@ -90,31 +91,55 @@ def bbox(self) -> List[float]:
             self._get_float("PROJECTION_ATTRIBUTES/CORNER_LL_LAT_PRODUCT"),
             self._get_float("PROJECTION_ATTRIBUTES/CORNER_LR_LAT_PRODUCT")
         ]
-
-        return [min(lons), min(lats), max(lons), max(lats)]
+        geod = Geod(ellps="WGS84")
+        offset = self.sr_gsd / 2
+        _, _, bottom_distance = geod.inv(lons[2], lats[2], lons[3], lats[3])
+        bottom_offset = offset * (lons[3] - lons[2]) / bottom_distance
+        _, _, top_distance = geod.inv(lons[0], lats[0], lons[1], lats[1])
+        top_offset = offset * (lons[1] - lons[0]) / top_distance
+        _, _, lat_distance = geod.inv(lons[0], lats[0], lons[2], lats[2])
+        lat_offset = offset * (lats[0] - lats[2]) / lat_distance
+        return [
+            min(lons) - bottom_offset,
+            min(lats) - lat_offset,
+            max(lons) + top_offset,
+            max(lats) + lat_offset
+        ]
 
     @property
     def proj_bbox(self) -> List[float]:
+        # USGS metadata provide bounds at the center of the pixel, but
+        # GDAL/rasterio transforms are to edge of pixel.
+        # https://github.com/stac-utils/stactools/issues/117
+        offset = self.sr_gsd / 2
         xs = [
             self._get_float(
-                "PROJECTION_ATTRIBUTES/CORNER_UL_PROJECTION_X_PRODUCT"),
+                "PROJECTION_ATTRIBUTES/CORNER_UL_PROJECTION_X_PRODUCT") -
+            offset,
             self._get_float(
-                "PROJECTION_ATTRIBUTES/CORNER_UR_PROJECTION_X_PRODUCT"),
+                "PROJECTION_ATTRIBUTES/CORNER_UR_PROJECTION_X_PRODUCT") +
+            offset,
             self._get_float(
-                "PROJECTION_ATTRIBUTES/CORNER_LL_PROJECTION_X_PRODUCT"),
+                "PROJECTION_ATTRIBUTES/CORNER_LL_PROJECTION_X_PRODUCT") -
+            offset,
             self._get_float(
-                "PROJECTION_ATTRIBUTES/CORNER_LR_PROJECTION_X_PRODUCT")
+                "PROJECTION_ATTRIBUTES/CORNER_LR_PROJECTION_X_PRODUCT") +
+            offset
         ]
 
         ys = [
             self._get_float(
-                "PROJECTION_ATTRIBUTES/CORNER_UL_PROJECTION_Y_PRODUCT"),
+                "PROJECTION_ATTRIBUTES/CORNER_UL_PROJECTION_Y_PRODUCT") +
+            offset,
             self._get_float(
-                "PROJECTION_ATTRIBUTES/CORNER_UR_PROJECTION_Y_PRODUCT"),
+                "PROJECTION_ATTRIBUTES/CORNER_UR_PROJECTION_Y_PRODUCT") +
+            offset,
             self._get_float(
-                "PROJECTION_ATTRIBUTES/CORNER_LL_PROJECTION_Y_PRODUCT"),
+                "PROJECTION_ATTRIBUTES/CORNER_LL_PROJECTION_Y_PRODUCT") -
+            offset,
             self._get_float(
-                "PROJECTION_ATTRIBUTES/CORNER_LR_PROJECTION_Y_PRODUCT")
+                "PROJECTION_ATTRIBUTES/CORNER_LR_PROJECTION_Y_PRODUCT") -
+            offset
         ]
 
         return [min(xs), min(ys), max(xs), max(ys)]

stactools_sentinel2/stactools/sentinel2/cog.py

@@ -2,11 +2,13 @@
 import os
 import shutil
 from tempfile import TemporaryDirectory
-from subprocess import Popen, PIPE, STDOUT
 
 import pystac
 from pystac.utils import make_absolute_href
 
+from stactools.core.utils.subprocess import call
+from stactools.core.utils.convert import cogify
+
 logger = logging.getLogger(__name__)
 
 # TODO: Refactor this into general-purpose cogification of item asset
@@ -52,23 +54,5 @@ def is_non_cog_image(asset):
         or asset.media_type == pystac.MediaType.JPEG2000)
 
 
-def call(command):
-    def log_subprocess_output(pipe):
-        for line in iter(pipe.readline, b''):  # b'\n'-separated lines
-            logger.info(line.decode("utf-8").strip('\n'))
-
-    process = Popen(command, stdout=PIPE, stderr=STDOUT)
-    with process.stdout:
-        log_subprocess_output(process.stdout)
-    return process.wait()  # 0 means success
-
-
-def cogify(input_path, output_path):
-    call([
-        'gdal_translate', '-of', 'COG', '-co', 'compress=deflate', input_path,
-        output_path
-    ])
-
-
 def reproject(input_path, output_path):
     call(['gdalwarp', '-t_srs', 'epsg:3857', input_path, output_path])

stactools_sentinel2/stactools/sentinel2/constants.py

@@ -93,3 +93,57 @@
                 center_wavelength=2.190,
                 full_width_half_max=0.242),
 }
+
+# A dict describing the resolutions that are
+# available for each band as separate assets.
+# The first resolution is the sensor gsd; others
+# are downscaled versions.
+BANDS_TO_RESOLUTIONS = {
+    'B01': [
+        60,
+    ],
+    'B02': [
+        10,
+        20,
+        60,
+    ],
+    'B03': [
+        10,
+        20,
+        60,
+    ],
+    'B04': [
+        10,
+        20,
+        60,
+    ],
+    'B05': [
+        20,
+        20,
+        60,
+    ],
+    'B06': [
+        20,
+        60,
+    ],
+    'B07': [
+        20,
+        60,
+    ],
+    'B08': [
+        10,
+    ],
+    'B8A': [
+        20,
+        60,
+    ],
+    'B09': [60],
+    'B11': [
+        20,
+        60,
+    ],
+    'B12': [
+        20,
+        60,
+    ],
+}

stactools_sentinel2/stactools/sentinel2/granule_metadata.py

@@ -149,7 +149,7 @@ def metadata_dict(self):
                 float,
                 self._image_content_node.find_text(
                     'HIGH_PROBA_CLOUDS_PERCENTAGE')),
-            's2:thin_sirrus_percentage':
+            's2:thin_cirrus_percentage':
             map_opt(
                 float,
                 self._image_content_node.find_text('THIN_CIRRUS_PERCENTAGE')),

stactools_sentinel2/stactools/sentinel2/stac.py

@@ -12,12 +12,10 @@
 from stactools.sentinel2.product_metadata import ProductMetadata
 from stactools.sentinel2.granule_metadata import GranuleMetadata
 from stactools.sentinel2.utils import extract_gsd
-from stactools.sentinel2.constants import (DATASTRIP_METADATA_ASSET_KEY,
-                                           SENTINEL_PROVIDER, SENTINEL_LICENSE,
-                                           SENTINEL_BANDS,
-                                           SENTINEL_INSTRUMENTS,
-                                           SENTINEL_CONSTELLATION,
-                                           INSPIRE_METADATA_ASSET_KEY)
+from stactools.sentinel2.constants import (
+    BANDS_TO_RESOLUTIONS, DATASTRIP_METADATA_ASSET_KEY, SENTINEL_PROVIDER,
+    SENTINEL_LICENSE, SENTINEL_BANDS, SENTINEL_INSTRUMENTS,
+    SENTINEL_CONSTELLATION, INSPIRE_METADATA_ASSET_KEY)
 
 logger = logging.getLogger(__name__)
 
@@ -172,29 +170,52 @@ def image_asset_from_href(
         return ('preview', asset)
 
     # Extract gsd and proj info
-    gsd = extract_gsd(asset_href)
-    shape = list(resolution_to_shape[int(gsd)])
+    filename_gsd = extract_gsd(asset_href)
+    shape = list(resolution_to_shape[int(filename_gsd)])
     transform = transform_from_bbox(proj_bbox, shape)
 
-    def set_asset_properties(asset):
-        item.common_metadata.set_gsd(gsd, asset)
+    def set_asset_properties(asset: pystac.Asset,
+                             band_gsd: Optional[int] = None):
+        if band_gsd:
+            item.common_metadata.set_gsd(band_gsd, asset)
         item.ext.projection.set_shape(shape, asset)
         item.ext.projection.set_bbox(proj_bbox, asset)
         item.ext.projection.set_transform(transform, asset)
 
     # Handle band image
 
-    band_id_search = re.search(r'_(B\w{2})_', asset_href)
+    band_id_search = re.search(r'_(B\w{2})', asset_href)
     if band_id_search is not None:
-        band_id = band_id_search.group(1)
+        band_id, href_res = os.path.splitext(asset_href)[0].split('_')[-2:]
         band = SENTINEL_BANDS[band_id]
+
+        # Get the asset resolution from the file name.
+        # If the asset resolution is the band GSD, then
+        # include the gsd information for that asset. Otherwise,
+        # do not include the GSD information in the asset
+        # as this may be confusing for users given that the
+        # raster spatial resolution and gsd will differ.
+        # See https://github.com/radiantearth/stac-spec/issues/1096
+        asset_res = int(href_res.replace('m', ''))
+        band_gsd: Optional[int] = None
+        if asset_res == BANDS_TO_RESOLUTIONS[band_id][0]:
+            asset_key = band_id
+            band_gsd = asset_res
+        else:
+            # If this isn't the default resolution, use the raster
+            # resolution in the asset key.
+            # TODO: Use the raster extension and spatial_resolution
+            # property to encode the spatial resolution of all assets.
+            asset_key = f'{band_id}_{asset_res}m'
+
         asset = pystac.Asset(href=asset_href,
                              media_type=asset_media_type,
-                             title=band.description,
+                             title=f'{band.description} - {href_res}',
                              roles=['data'])
+
         item.ext.eo.set_bands([SENTINEL_BANDS[band_id]], asset)
-        set_asset_properties(asset)
-        return (band_id, asset)
+        set_asset_properties(asset, band_gsd=band_gsd)
+        return (asset_key, asset)
 
     # Handle auxiliary images