Merge pull request #95 from ACCESS-Cloud-Based-InSAR/dev

v2.5.7

Author: cmarshak

.github/workflows/changelog.yml

@@ -13,4 +13,4 @@ on:
 
 jobs:
   call-changelog-check-workflow:
-    uses: ASFHyP3/actions/.github/workflows/[email protected].1
+    uses: ASFHyP3/actions/.github/workflows/[email protected].2

.github/workflows/labeled-pr.yml

@@ -12,4 +12,4 @@ on:
 
 jobs:
   call-labeled-pr-check-workflow:
-    uses: ASFHyP3/actions/.github/workflows/[email protected].1
+    uses: ASFHyP3/actions/.github/workflows/[email protected].2

.github/workflows/release.yml

@@ -7,7 +7,7 @@ on:
 
 jobs:
   call-release-workflow:
-    uses: ASFHyP3/actions/.github/workflows/[email protected].1
+    uses: ASFHyP3/actions/.github/workflows/[email protected].2
     with:
       release_prefix: dem-stitcher
       develop_branch: dev

.github/workflows/static-analysis.yml

@@ -4,12 +4,12 @@ on: push
 
 jobs:
   call-flake8-workflow:
-    uses: ASFHyP3/actions/.github/workflows/[email protected].1
+    uses: ASFHyP3/actions/.github/workflows/[email protected].2
     with:
       local_package_names: dem_stitcher
 
   call-secrets-analysis-workflow:
-    uses: ASFHyP3/actions/.github/workflows/[email protected].1
+    uses: ASFHyP3/actions/.github/workflows/[email protected].2
 
   call-ruff-workflow:
-    uses: ASFHyP3/actions/.github/workflows/[email protected].1
+    uses: ASFHyP3/actions/.github/workflows/[email protected].2

.github/workflows/tag.yml

@@ -7,7 +7,7 @@ on:
 
 jobs:
   call-bump-version-workflow:
-    uses: ASFHyP3/actions/.github/workflows/[email protected].1
+    uses: ASFHyP3/actions/.github/workflows/[email protected].2
     with:
       user: access-cloud-insar-team
       email: [email protected]

CHANGELOG.md

@@ -6,6 +6,16 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/)
 and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [2.5.7]
+
+### Fixed
+* Check for Earthdata credentials in netrc (adapted from Joe Kennedy/Forrest Williams) resolving isse [#83](https://github.com/ACCESS-Cloud-Based-InSAR/dem-stitcher/issues/83)
+  * when no credentials in netrc are present when requesting data for `nasadem` or `srtm_v3`, there is a human readable error instructing user to update their `~/.netrc`.
+* Updates some ruff linting
+  * Ensures ruff in `environment.yml`
+  * Ensure single quotes for consistency. 
+
+
 ## [2.5.6]
 * Updated URLs for downloading geoids from agisoft.com. Fixes [#88](https://github.com/ACCESS-Cloud-Based-InSAR/dem-stitcher/issues/88).
 

README.md

@@ -67,8 +67,7 @@ The creation metadata unrelated to georeferencing (e.g. the `compress` key or va
 ## Credentials
 
 The accessing of NASADEM and SRTM require earthdata login credentials to be put into the `~/.netrc` file. If these are not present, the stitcher will
-fail with `BadZipFile Error` as we use `requests` to obtain zipped data and load the data using `rasterio`. An entry in the `.netrc` will look like:
-
+fail with `ValueError` asking you to update the `~/.netrc`. The appropriate entry appears as:
 ```
 machine urs.earthdata.nasa.gov
     login <username>

dem_stitcher/__init__.py

@@ -11,17 +11,17 @@
 except PackageNotFoundError:
     __version__ = None
     warnings.warn(
-        "package is not installed!\n"
-        "Install in editable/develop mode via (from the top of this repo):\n"
-        "   python -m pip install -e .\n",
+        'package is not installed!\n'
+        'Install in editable/develop mode via (from the top of this repo):\n'
+        '   python -m pip install -e .\n',
         RuntimeWarning,
     )
 
 
 __all__ = [
-    "get_dem_tile_paths",
-    "get_global_dem_tile_extents",
-    "get_overlapping_dem_tiles",
-    "stitch_dem",
-    "__version__",
+    'get_dem_tile_paths',
+    'get_global_dem_tile_extents',
+    'get_overlapping_dem_tiles',
+    'stitch_dem',
+    '__version__',
 ]

dem_stitcher/credentials.py

@@ -0,0 +1,23 @@
+import netrc
+from pathlib import Path
+
+
+def ensure_earthdata_credentials(
+    host: str = 'urs.earthdata.nasa.gov',
+):
+    """
+    Source: DockerizedTopsapp / Authors: Joseph Kennedy, Forrest Williams, and Andrew Johnston
+
+    Ensures Earthdata credentials are provided in ~/.netrc
+
+    Earthdata username and password may be provided by, in order of preference, one of:
+       * `netrc_file`
+       * `username` and `password`
+    and will be written to the ~/.netrc file if it doesn't already exist.
+    """
+    netrc_file = Path.home() / '.netrc'
+    try:
+        dot_netrc = netrc.netrc(netrc_file)
+        _, _, _ = dot_netrc.authenticators(host)
+    except (FileNotFoundError, netrc.NetrcParseError, TypeError):
+        raise ValueError(f'Please provide valid Earthdata login credentials via {netrc_file}')

dem_stitcher/datasets.py

@@ -12,11 +12,11 @@
 from .exceptions import DEMNotSupported
 from .geojson_io import read_geojson_gzip
 
-DATA_PATH = Path(__file__).parents[0].absolute() / "data"
+DATA_PATH = Path(__file__).parents[0].absolute() / 'data'
 
 # Get Datasets
-_DATASET_PATHS = list(DATA_PATH.glob("*.geojson.zip"))
-DATASETS = list(map(lambda x: x.name.split(".")[0], _DATASET_PATHS))
+_DATASET_PATHS = list(DATA_PATH.glob('*.geojson.zip'))
+DATASETS = list(map(lambda x: x.name.split('.')[0], _DATASET_PATHS))
 
 
 def get_available_datasets():
@@ -45,8 +45,8 @@ def get_global_dem_tile_extents(dataset: str) -> gpd.GeoDataFrame:
     """
     if dataset not in DATASETS:
         raise DEMNotSupported(f'{dataset} must be in {", ".join(DATASETS)}')
-    df = read_geojson_gzip(DATA_PATH / f"{dataset}.geojson.zip")
-    df["dem_name"] = dataset
+    df = read_geojson_gzip(DATA_PATH / f'{dataset}.geojson.zip')
+    df['dem_name'] = dataset
     df.crs = CRS.from_epsg(4326)
     return df
 
@@ -81,9 +81,9 @@ def get_overlapping_dem_tiles(bounds: list, dem_name: str) -> gpd.GeoDataFrame:
     crossing = get_dateline_crossing(bounds)
     if crossing:
         warn(
-            "Getting tiles across dateline on the opposite hemisphere; "
-            f"The source tiles will be {- 2 * crossing} deg along the"
-            "longitudinal axis from the extent requested",
+            'Getting tiles across dateline on the opposite hemisphere; '
+            f'The source tiles will be {- 2 * crossing} deg along the'
+            'longitudinal axis from the extent requested',
             category=UserWarning,
         )
         df_tiles_all_translated = df_tiles_all.copy()
@@ -100,18 +100,18 @@ def get_overlapping_dem_tiles(bounds: list, dem_name: str) -> gpd.GeoDataFrame:
     if not df_tiles.empty:
         # Degenerate geometries raise warning in shapely - intersection is black box to us
         with warnings.catch_warnings():
-            warnings.simplefilter("ignore", category=RuntimeWarning)
+            warnings.simplefilter('ignore', category=RuntimeWarning)
             df_tiles_intersection = df_tiles.geometry.intersection(box_geo)
-            geo_type_index = df_tiles_intersection.geometry.map(lambda geo: geo.geom_type == "Polygon")
+            geo_type_index = df_tiles_intersection.geometry.map(lambda geo: geo.geom_type == 'Polygon')
         df_tiles = df_tiles[geo_type_index].copy()
 
     # Merging is order dependent - ensures consistency
-    df_tiles = df_tiles.sort_values(by="tile_id")
+    df_tiles = df_tiles.sort_values(by='tile_id')
     df_tiles = df_tiles.reset_index(drop=True)
     return df_tiles
 
 
 def intersects_missing_glo_30_tiles(extent: list) -> bool:
     extent_geo = box(*extent)
-    df_missing = get_overlapping_dem_tiles(extent, "glo_90_missing")
+    df_missing = get_overlapping_dem_tiles(extent, 'glo_90_missing')
     return df_missing.intersects(extent_geo).sum() > 0

dem_stitcher/dateline.py

@@ -10,19 +10,19 @@ def check_4326_bounds(bounds: list) -> bool:
     xmin, ymin, xmax, ymax = bounds
 
     if (xmin > xmax) or (ymin > ymax):
-        raise Incorrect4326Bounds("Ensure xmin <= xmax and ymin <= ymax")
+        raise Incorrect4326Bounds('Ensure xmin <= xmax and ymin <= ymax')
 
     standard_4326_box = box(-180, -90, 180, 90)
     bounds_box = box(*bounds)
 
     if not (standard_4326_box.intersects(bounds_box)):
         raise Incorrect4326Bounds(
-            "Make sure bounds have intersection over standard 4326 CRS i.e. "
-            "between longitude -180 and 180 and latitude -90 and 90."
+            'Make sure bounds have intersection over standard 4326 CRS i.e. '
+            'between longitude -180 and 180 and latitude -90 and 90.'
         )
 
     if (ymin < -90) or (ymax > 90):
-        raise Incorrect4326Bounds("Boxes beyond the North/South Pole at +/- 90 Latitude not supported")
+        raise Incorrect4326Bounds('Boxes beyond the North/South Pole at +/- 90 Latitude not supported')
 
     return True
 
@@ -65,7 +65,7 @@ def get_dateline_crossing(bounds: list) -> int:
         return 180
 
     elif (xmin <= -180) and (xmax >= 180):
-        raise DoubleDatelineCrossing("Shrink your bounding area")
+        raise DoubleDatelineCrossing('Shrink your bounding area')
 
 
 def split_extent_across_dateline(extent: list) -> tuple[list]:
@@ -99,7 +99,7 @@ def split_extent_across_dateline(extent: list) -> tuple[list]:
         multipolygon = extent_box.union(extent_box_t)
 
         with warnings.catch_warnings():
-            warnings.simplefilter("ignore", category=RuntimeWarning)
+            warnings.simplefilter('ignore', category=RuntimeWarning)
             bounds_l = list(multipolygon.intersection(left_hemisphere).bounds)
             bounds_r = list(multipolygon.intersection(right_hemisphere).bounds)
         return (bounds_l, bounds_r)

dem_stitcher/dem_readers.py

@@ -15,9 +15,9 @@ def read_dem(dem_path: str) -> rasterio.DatasetReader:
     return dem_arr, dem_profile
 
 
-def read_dem_bytes(dem_path: str, suffix: str = ".img") -> bytes:
+def read_dem_bytes(dem_path: str, suffix: str = '.img') -> bytes:
     # online
-    if (dem_path[:7] == "http://") or (dem_path[:8] == "https://"):
+    if (dem_path[:7] == 'http://') or (dem_path[:8] == 'https://'):
         resp = requests.get(dem_path)
         data = io.BytesIO(resp.content)
     # local file
@@ -35,17 +35,17 @@ def read_dem_bytes(dem_path: str, suffix: str = ".img") -> bytes:
     return img_bytes
 
 
-def read_srtm(dem_path: str, version="srtm") -> Tuple[np.ndarray, dict]:
-    img_bytes = read_dem_bytes(dem_path, suffix=".hgt")
+def read_srtm(dem_path: str, version='srtm') -> Tuple[np.ndarray, dict]:
+    img_bytes = read_dem_bytes(dem_path, suffix='.hgt')
     # The gdal driver hgt depends on filename convention
-    filename = dem_path.split("/")[-1]
-    if version == "srtm":
-        filename = filename.replace(".zip", "")
-    elif version == "nasadem":
-        filename = filename.replace(".zip", ".hgt")
-        filename = filename.replace("NASADEM_HGT_", "")
+    filename = dem_path.split('/')[-1]
+    if version == 'srtm':
+        filename = filename.replace('.zip', '')
+    elif version == 'nasadem':
+        filename = filename.replace('.zip', '.hgt')
+        filename = filename.replace('NASADEM_HGT_', '')
     else:
-        raise ValueError("version must be either nasadem or srtm")
+        raise ValueError('version must be either nasadem or srtm')
 
     with MemoryFile(img_bytes, filename=filename) as memfile:
         with memfile.open() as dataset:
@@ -56,4 +56,4 @@ def read_srtm(dem_path: str, version="srtm") -> Tuple[np.ndarray, dict]:
 
 
 def read_nasadem(dem_path: str) -> Tuple[np.ndarray, dict]:
-    return read_srtm(dem_path, version="nasadem")
+    return read_srtm(dem_path, version='nasadem')

dem_stitcher/geoid.py

@@ -11,11 +11,11 @@
 from .rio_tools import reproject_arr_to_match_profile, translate_profile
 from .rio_window import read_raster_from_window
 
-AGISOFT_URL = "https://s3-eu-west-1.amazonaws.com/download.agisoft.com/geoids"
+AGISOFT_URL = 'https://s3-eu-west-1.amazonaws.com/download.agisoft.com/geoids'
 GEOID_PATHS_AGI = {
-    "geoid_18": f"{DATA_PATH}/geoid_18.tif",
-    "egm_08": f"{AGISOFT_URL}/egm2008-1.tif",
-    "egm_96": f"{AGISOFT_URL}/egm96-15.tif",
+    'geoid_18': f'{DATA_PATH}/geoid_18.tif',
+    'egm_08': f'{AGISOFT_URL}/egm2008-1.tif',
+    'egm_96': f'{AGISOFT_URL}/egm96-15.tif',
 }
 
 
@@ -46,7 +46,7 @@ def read_geoid(geoid_name: str, extent: list = None, res_buffer: int = 1) -> tup
             geoid_arr_r, geoid_profile_r = read_raster_from_window(
                 geoid_path, extent_r, extent_crs, res_buffer=res_buffer
             )
-            res_x = geoid_profile_l["transform"].a
+            res_x = geoid_profile_l['transform'].a
             if crossing == 180:
                 geoid_profile_l = translate_profile(geoid_profile_l, 360 / res_x, 0)
             else:
@@ -55,42 +55,42 @@ def read_geoid(geoid_name: str, extent: list = None, res_buffer: int = 1) -> tup
                 [geoid_arr_l, geoid_arr_r], [geoid_profile_l, geoid_profile_r]
             )
     # Transform nodata to nan
-    geoid_arr = geoid_arr.astype("float32")
-    geoid_arr[geoid_profile["nodata"] == geoid_arr] = np.nan
-    geoid_profile["nodata"] = np.nan
+    geoid_arr = geoid_arr.astype('float32')
+    geoid_arr[geoid_profile['nodata'] == geoid_arr] = np.nan
+    geoid_profile['nodata'] = np.nan
 
     return geoid_arr, geoid_profile
 
 
 def remove_geoid(
-    dem_arr: np.ndarray, dem_profile: dict, geoid_name: str, dem_area_or_point: str = "Area", res_buffer: int = 2
+    dem_arr: np.ndarray, dem_profile: dict, geoid_name: str, dem_area_or_point: str = 'Area', res_buffer: int = 2
 ) -> np.ndarray:
-    assert dem_area_or_point in ["Point", "Area"]
+    assert dem_area_or_point in ['Point', 'Area']
 
-    extent = array_bounds(dem_profile["height"], dem_profile["width"], dem_profile["transform"])
+    extent = array_bounds(dem_profile['height'], dem_profile['width'], dem_profile['transform'])
 
     geoid_arr, geoid_profile = read_geoid(geoid_name, extent=list(extent), res_buffer=res_buffer)
 
-    t_dem = dem_profile["transform"]
-    t_geoid = geoid_profile["transform"]
+    t_dem = dem_profile['transform']
+    t_geoid = geoid_profile['transform']
     res_dem = max(t_dem.a, abs(t_dem.e))
     res_geoid = max(t_geoid.a, abs(t_geoid.e))
 
     if res_geoid * res_buffer <= res_dem:
         buffer_recommendation = int(np.ceil(res_dem / res_geoid))
         warning = (
-            "The dem resolution is larger than the geoid resolution and its buffer; "
-            "Edges resampled with bilinear interpolation will be inconsistent so select larger buffer."
-            f"Select a `res_buffer = {buffer_recommendation}`"
+            'The dem resolution is larger than the geoid resolution and its buffer; '
+            'Edges resampled with bilinear interpolation will be inconsistent so select larger buffer.'
+            f'Select a `res_buffer = {buffer_recommendation}`'
         )
         warnings.warn(warning, category=UserWarning)
 
     # Translate geoid if necessary as all geoids have Area tag
-    if dem_area_or_point == "Point":
+    if dem_area_or_point == 'Point':
         shift = -0.5
         geoid_profile = translate_profile(geoid_profile, shift, shift)
 
-    geoid_offset, _ = reproject_arr_to_match_profile(geoid_arr, geoid_profile, dem_profile, resampling="bilinear")
+    geoid_offset, _ = reproject_arr_to_match_profile(geoid_arr, geoid_profile, dem_profile, resampling='bilinear')
 
     dem_arr_offset = dem_arr + geoid_offset
     return dem_arr_offset

dem_stitcher/geojson_io.py

@@ -9,26 +9,26 @@
 
 
 def read_geojson_gzip(input_zip_path: Union[str, Path]) -> gpd.GeoDataFrame:
-    with gzip.GzipFile(input_zip_path, "r") as file_in:
-        data_gjson = json.loads(file_in.read().decode("utf-8"))
-    return gpd.GeoDataFrame.from_features(data_gjson["features"], crs=CRS.from_epsg(4326))
+    with gzip.GzipFile(input_zip_path, 'r') as file_in:
+        data_gjson = json.loads(file_in.read().decode('utf-8'))
+    return gpd.GeoDataFrame.from_features(data_gjson['features'], crs=CRS.from_epsg(4326))
 
 
 def to_geojson_obj(geodataframe: gpd.geodataframe.GeoDataFrame) -> dict:
-    features = geodataframe.to_dict("records")
+    features = geodataframe.to_dict('records')
 
     def mapping_geojson(entry):
-        geometry = entry.pop("geometry")
-        new_entry = {"type": "Feature", "properties": entry, "geometry": shapely.geometry.mapping(geometry)}
+        geometry = entry.pop('geometry')
+        new_entry = {'type': 'Feature', 'properties': entry, 'geometry': shapely.geometry.mapping(geometry)}
         return new_entry
 
     features = list(map(mapping_geojson, features))
-    geojson = {"type": "FeatureCollection", "features": features}
+    geojson = {'type': 'FeatureCollection', 'features': features}
     return geojson
 
 
 def to_geojson_gzip(geodataframe: gpd.geodataframe.GeoDataFrame, dest_path: str) -> Path:
     geojson_ob = to_geojson_obj(geodataframe)
-    with gzip.GzipFile(dest_path, "w") as file_out:
-        file_out.write(json.dumps(geojson_ob).encode("utf-8"))
+    with gzip.GzipFile(dest_path, 'w') as file_out:
+        file_out.write(json.dumps(geojson_ob).encode('utf-8'))
     return dest_path