Add support for USGS 3DEP (formerly NED) items and collections

.codespellignore

@@ -0,0 +1 @@
+ned

.readthedocs.yml

@@ -31,6 +31,8 @@ python:
       path: stactools_landsat/
     - method: pip
       path: stactools_planet/
+    - method: pip
+      path: stactools_threedep/
     - method: pip
       path: stactools_browse/
     - method: pip

README.md

@@ -79,6 +79,7 @@ See the [documentation page](https://stactools.readthedocs.io/en/latest/) for th
 | `stactools_planet`            | Methods and commands for working with planet data                |
 | `stactools_landsat`           | Methods and commands for working with landsat data        |
 | `stactools_sentinel2`           | Methods and commands for working with Sentinel 2 data        |
+| `stactools_threedep`          | Methods and commands for working with 3DEP (formerly NED) elevation data        |
 | `stactools_browse`            | Contains a command for launching stac-browser against a local STAC |
 
 Subpackages are symlinked to the `stactools` directory in this repo to allow them to be importable for python running at the top level directory of the repository clone.

docker/Dockerfile

@@ -82,6 +82,10 @@ RUN pip install -r /tmp/planet/requirements.txt
 COPY stactools_sentinel2/requirements.txt /tmp/sentinel2/requirements.txt
 RUN pip install -r /tmp/sentinel2/requirements.txt
 
+# 3DEP
+COPY stactools_threedep/requirements.txt /tmp/threedep/requirements.txt
+RUN pip install -r /tmp/threedep/requirements.txt
+
 # Jupyter
 RUN pip install jupyter==1.0.0
 

scripts/env

@@ -9,6 +9,7 @@ STACTOOLS_SUBPACKAGE_DIRS=(
     "stactools_landsat"
     "stactools_planet"
     "stactools_naip"
+    "stactools_threedep"
     "stactools_browse"
     "stactools_sentinel2"
 );
@@ -22,6 +23,7 @@ STACTOOLS_COVERAGE_DIRS=(
     "stactools_landsat/stactools/landsat"
     "stactools_planet/stactools/planet"
     "stactools_naip/stactools/naip"
+    "stactools_threedep/stactools/threedep"
     "stactools_browse/stactools/browse"
     "stactools_sentinel2/stactools/sentinel2"
 );

setup.py

@@ -34,6 +34,7 @@ def requirement_name(self):
     Subpackage('landsat', is_extra=True),
     Subpackage('naip', is_extra=True),
     Subpackage('planet', is_extra=True),
+    Subpackage('threedep', is_extra=True),
     Subpackage('browse', is_extra=True),
     Subpackage('sentinel2', is_extra=True)
 ]

stactools/threedep

@@ -0,0 +1 @@
+../stactools_threedep/stactools/threedep/

stactools_threedep/README.md

@@ -0,0 +1,46 @@
+# stactools.threedep
+
+The 3D Elevation Program (3DEP), formerly known as the National Elevation Dataset (NED), is elevation data for the United States.
+Because `3dep` isn't a valid Python package name, this package is named `stactools.threedep`.
+
+## Usage
+
+3DEP items are identified by two attributes: `product` and `id`.
+`product` is a string that corresponds to the nominal resolution of the data.
+Options are:
+
+- "1": 1 arc-second DEMs
+- "13" 1/3 arc-second DEMs
+
+`id` is a lat-lon identifier, e.g. `n41w106`.
+
+### Command line
+
+To build a `pystac.Collection` directly from AWS into a directory named `usgs-3dep-stac`:
+
+```bash
+stac threedep create-collection usgs-3dep-stac
+```
+
+If you want to two-step the process, you can download the metadata first:
+
+```bash
+stac threedep download-metadata usgs-3dep-metadata
+stac threedep create-collection usgs-3dep-stac --source usgs-3dep-metadata
+```
+
+### API
+
+To create an item from AWS, use `stactools.threedep.stac.create_item`:
+
+```python
+from stactools.threedep import stac
+item = stac.create_item_from_product_and_id("1", "n41w106")
+```
+
+You can also create an item directly from the href of a metadata XML file anyhwere:
+
+```python
+from stactools.threedep import stac
+item = stac.create_item("ftp://rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/Elevation/1/TIFF/n41w106/USGS_1_n41w106.xml")
+```

stactools_threedep/requirements.txt

@@ -0,0 +1 @@
+boto3==1.17.40

stactools_threedep/setup.py

@@ -0,0 +1,47 @@
+import os
+from os.path import (basename, splitext)
+from imp import load_source
+from setuptools import setup, find_namespace_packages
+from glob import glob
+import io
+
+name = 'stactools_threedep'
+description = ("Subpackage for working with 3DEP data in stactools, "
+               "a command line tool and Python library for working with STAC.")
+
+__version__ = load_source(
+    'stactools.threedep.version',
+    os.path.join(os.path.dirname(__file__),
+                 'stactools/threedep/version.py')).__version__
+
+here = os.path.abspath(os.path.dirname(__file__))
+
+# get the dependencies and installs
+with io.open(os.path.join(here, 'requirements.txt'), encoding='utf-8') as f:
+    install_requires = [line.split(' ')[0] for line in f.read().split('\n')]
+
+# Add stactools subpackage dependencies
+install_requires.extend(['stactools_core=={}'.format(__version__)])
+
+with open(os.path.join(here, 'README.md')) as readme_file:
+    readme = readme_file.read()
+
+setup(name=name,
+      description=description,
+      version=__version__,
+      long_description=readme,
+      long_description_content_type="text/markdown",
+      author="Pete Gadomski",
+      author_email='[email protected]',
+      url='https://github.com/stac-utils/stactools.git',
+      packages=find_namespace_packages(),
+      py_modules=[
+          splitext(basename(path))[0] for path in glob('stactools/*.py')
+      ],
+      include_package_data=False,
+      install_requires=install_requires,
+      license="Apache Software License 2.0",
+      keywords=[
+          'stactools', 'psytac', '3DEP', 'NED', 'DEM', 'elevation', 'raster',
+          'catalog', 'STAC'
+      ])

stactools_threedep/stactools/threedep/__init__.py

@@ -0,0 +1,13 @@
+import stactools.core
+
+from stactools.threedep.metadata import Metadata
+
+stactools.core.use_fsspec()
+
+
+def register_plugin(registry):
+    from stactools.threedep import commands
+    registry.register_subcommand(commands.create_threedep_command)
+
+
+__all__ = [Metadata]

stactools_threedep/stactools/threedep/commands.py

@@ -0,0 +1,126 @@
+import click
+import os.path
+
+from pystac import Catalog, Collection, Extent, CatalogType, STAC_IO
+
+from stactools.threedep.constants import (PRODUCTS, DESCRIPTION, USGS_PROVIDER,
+                                          USGS_FTP_BASE, USGS_3DEP_ID)
+from stactools.threedep import utils, stac
+
+
+def create_threedep_command(cli):
+    """Creates the threedep command line utility."""
+    @cli.group("threedep", short_help="Work with USGS 3DEP elevation data.")
+    def threedep():
+        pass
+
+    @threedep.command(
+        "create-collection",
+        short_help="Create a STAC collection for existing USGS 3DEP data")
+    @click.argument("destination")
+    @click.option("-s",
+                  "--source",
+                  help="The href of a directory tree containing metadata",
+                  default=None)
+    @click.option("-i",
+                  "--id",
+                  multiple=True,
+                  help="Ids to fetch. If not provided, will fetch all IDs.")
+    @click.option("--quiet/--no-quiet", default=False)
+    def create_collection_command(destination, source, id, quiet):
+        """Creates a 3DEP collection in DESTINATION.
+
+        If SOURCE is not provided, will use the metadata in AWS. SOURCE is
+        expected to be a directory tree mirroring the structure on USGS, so
+        it is best created using `stac threedep download-metadata`.
+        """
+        base_ids = id  # not sure how to rename arguments in click
+        items = {}
+        for product in PRODUCTS:
+            items[product] = []
+            if base_ids:
+                ids = base_ids
+            else:
+                ids = utils.fetch_ids(product)
+            for id in ids:
+                item = stac.create_item_from_product_and_id(
+                    product, id, source)
+                items[product].append(item)
+                if not quiet:
+                    print(item.id)
+        all_items = [item for sublist in items.values() for item in sublist]
+        if not all_items:
+            raise Exception("no items found")
+        extent = Extent.from_items(all_items)
+        collection = Collection(
+            id=USGS_3DEP_ID,
+            title="USGS 3DEP DEMs",
+            keywords=["USGS", "3DEP", "NED", "DEM", "elevation"],
+            providers=[USGS_PROVIDER],
+            description=DESCRIPTION,
+            extent=extent,
+            license="PDDL-1.0")
+        for product in PRODUCTS:
+            if product == "1":
+                title = "1 arc-second"
+                description = "USGS 3DEP 1 arc-second DEMs"
+            elif product == "13":
+                title = "1/3 arc-second"
+                description = "USGS 3DEP 1/3 arc-second DEMs"
+            else:
+                raise NotImplementedError
+            catalog = Catalog(id=product, description=description, title=title)
+            collection.add_child(catalog)
+            catalog.add_items(items[product])
+        collection.update_extent_from_items()
+        collection.normalize_hrefs(destination)
+        collection.save(catalog_type=CatalogType.SELF_CONTAINED)
+        collection.validate()
+
+    @threedep.command("download-metadata",
+                      short_help="Download all metadata for USGS 3DEP data")
+    @click.argument("destination")
+    @click.option("-i",
+                  "--id",
+                  multiple=True,
+                  help="Ids to fetch. If not provided, will fetch all IDs.")
+    @click.option("--quiet/--no-quiet", default=False)
+    def download_metadata_command(destination, id, quiet):
+        """Creates a 3DEP collection in DESTINATION."""
+        base_ids = id  # not sure how to rename arguments in click
+        for product in PRODUCTS:
+            if base_ids:
+                ids = base_ids
+            else:
+                ids = utils.fetch_ids(product)
+            for id in ids:
+                path = utils.path(product,
+                                  id,
+                                  extension="xml",
+                                  base=destination)
+                if os.path.exists(path):
+                    if not quiet:
+                        print("{} exists, skipping download...".format(path))
+                    continue
+                os.makedirs(os.path.dirname(path), exist_ok=True)
+                source_path = utils.path(product,
+                                         id,
+                                         extension="xml",
+                                         base=USGS_FTP_BASE)
+                if not quiet:
+                    print("{} -> {}".format(source_path, path))
+                text = STAC_IO.read_text(source_path)
+                with open(path, "w") as f:
+                    f.write(text)
+
+    @threedep.command(
+        "fetch-ids",
+        short_help="Fetch all product ids and print them to stdout")
+    @click.argument("product")
+    @click.option("--usgs-ftp/--no-usgs-ftp",
+                  default=False,
+                  help="Fetch from the USGS FTP instead of AWS.")
+    def fetch_ids_command(product: str, usgs_ftp: bool):
+        """Fetches product ids and prints them to stdout."""
+        for id in utils.fetch_ids(product, use_usgs_ftp=usgs_ftp):
+            print(id)

stactools_threedep/stactools/threedep/constants.py

@@ -0,0 +1,28 @@
+# flake8: noqa
+
+from pyproj import CRS
+from pystac import Provider
+
+USGS_3DEP_ID = "usgs-3dep"
+THREEDEP_EPSG = 5498
+THREEDEP_CRS = CRS.from_epsg(THREEDEP_EPSG)
+PRODUCTS = ["1", "13"]
+LICENSE = "PDDL-1.0"
+
+DESCRIPTION = """The USGS 3D Elevation Program (3DEP) Datasets from The National Map are the primary elevation data product produced and distributed by the USGS. The 3DEP program provides a variety of resolution raster elevation data of the conterminous United States, Alaska, Hawaii, and the island territories. Some of the data sets such as the 1/3rd arc-second and 1 arc-second data set are derived from diverse source data sets that are processed to a specification with a consistent resolution, coordinate system, elevation units, and horizontal and vertical datums.  These seamless DEMs were referred to as the National Elevation Dataset (NED) from about 2000 through 2015 at which time they became the seamless DEM layers under the 3DEP program and the NED name and system were retired. Other 3DEP products include one-meter DEMs produced exclusively from high resolution light detection and ranging (lidar) source data and five-meter DEMs in Alaska as well as various source datasets including the lidar point cloud and interferometric synthetic aperture radar (Ifsar) digital surface models and intensity images. All 3DEP products are public domain. The 3DEP program is the logical result of the maturation of the long-standing USGS elevation program, which for many years concentrated on production of topographic map quadrangle-based digital elevation models. The 3DEP data  serves as the elevation layer of The National Map, and provides basic elevation information for earth science studies and mapping applications in the United States.
+
+The seamless DEM layers under the 3DEP program are a multi-resolution dataset that is updated continuously to integrate newly available, improved elevation source data. Seamless DEM data layers under the 3DEP program  are available nationally at grid spacings of 1 arc-second (approximately 30 meters) for the conterminous United States, and at 1/3, 1/9 arc-seconds  (approximately 10 and 3 meters, respectively) and 1 meter for parts of the United States. Most seamless DEM data for Alaska is available at 2-arc-second (about 60 meters) grid spacing, where only lower resolution source data exist. Part of Alaska is available at the 1/3 arc-second, 1 arc-second and 5 meter resolution. Efforts are continuing to have full coverage of the 5 meter elevation data over Alaska in the next couple of years.
+"""
+
+USGS_PROVIDER = Provider(
+    name="USGS",
+    roles=["producer", "processor", "host"],
+    url="https://www.usgs.gov/core-science-systems/ngp/3dep")
+
+USGS_FTP_SERVER = "rockyftp.cr.usgs.gov"
+USGS_FTP_BASE = f"ftp://{USGS_FTP_SERVER}/vdelivery/Datasets/Staged/Elevation"
+AWS_BUCKET = "prd-tnm"
+AWS_PREFIX = "StagedProducts/Elevation"
+AWS_BASE = f"https://{AWS_BUCKET}.s3.amazonaws.com/{AWS_PREFIX}"
+
+DEFAULT_BASE = AWS_BASE