access-cloud-based-insar/dem-stitcher

Download and merge DEM tiles

dem-stitcher's Issues

pip install for Python <=3.7 is v. 0.0.1

Platform:
Mac OS Big Sur, 11.6.4 (20G417) AND google Colab
Python v. 3.7

"pip install dem-stitcher"
Output:
Successfully installed dem-stitcher-0.0.1

It seems only 0.0.1 gets installed with Python <=3.7. Can you make this and/or a prior version compatible with at least Python=3? This is quite useful for Google Colab in particular, which only runs python 3.7.

Add Python 3.10 Back

Background

Some issues with gdal version 3.4 were noted here. Appears 3.5+ solves those issues.

Describe the solution you'd like

Pytest passes

Alternatives

N/A

Additional context

Want to have 3.10 via conda-forge

Pixel shift when dst_area_or_point == 'Point'

Hi,

Here is a slice of code in the stitcher.py:

    if dst_area_or_point == 'Point' and src_area_or_point == 'Area':
        x_shift = 1
        y_shift = 1
        array_shifted, profile_shifted = gdal_translate_profile(filepath, x_shift, y_shift)
    elif (dst_area_or_point == 'Area') and (src_area_or_point == 'Point'):
        shift = .5
        array_shifted, profile_shifted = gdal_translate_profile(filepath, shift, shift)
    # half shift down if glo30
    elif (dst_area_or_point == 'Point') and (src_area_or_point == 'Point'):
        x_shift = 1
        y_shift = 1
        array_shifted, profile_shifted = gdal_translate_profile(filepath, x_shift, y_shift)

I notice that x_shift and y_shift always equal to 1 if dst_area_or_point == 'Point' and no matter what src_area_or_point is. This looks a little strange to me. Is it a typo? If there is some reason behind it, I will appreciate it if you let me know.

Thanks!

Resampling and Profile Translations are Non-commutative and related Errors

The bug

Simple demonstration of bug:

If we specify extents contained entirely within a DEM tile and require no transformations (e.g. do not remove geoid and keep pixel center or area center coordinate system of the original tiles), we want dem-stitcher to return a subset of that tile that contains the extent. This is not the case. We will post a notebook to demonstrate this issue is resolved as well as an integration test.

In original version, we (well, the error was mine) used rasterio.merge (link) with the bounds keyword argument. I believed this function to be obtaining a window around the bound extents; however, as the rasterio source code reveals, there is resampling that is done to ensure that the extents fit nicely into the CRS, which may be beneficial in many, if not most, GIS applications.

More precisely, to go from Pixel center to Area center coordinates, we need to translate the original rasters and then perform any necessary resampling. If resampling is done before a translation of the geo-transform, then things go wrong (majorly).

Further, the window operations were also resampling under the hood. All these issues need to be fixed.

In a general sense, if we have two affine transformations T_r and T_t representing "resampling" and "translation" resepectively, then T_r * T_t is not the same as T_t * T_r (it's worth checking this on simple examples i.e. seeing that both the geo-transform/affine transformaion and the resulting rasters are different).

This was further complicated when we wanted to change the resolution at higher latitudes because glo-30 shrinks latitude spacing as documented here. Thus, more resampling was done.

Which order is correct? We will take the approach T_r * T_t i.e.:

1. First translating profiles (if required) on raw data
2. Resampling (if required) at the end of all transformations.

To Reproduce

To be filled in later.

Additional context

The issue #31 points out is due the fact this resampling is done so a shift was selected based on investigating sites. However, this is not a totally correct approach.

Work to be done

We need tests to make sure that the above bug is removed.

Install via Conda-Forge

Background

Install via conda-forge to ensure environment restrictions and more reliable for DockerizedTopsApp.

Describe the solution you'd like

In readme, conda install -c conda-forge dem_stitcher

Alternatives

Not applicable

Additional context

environment.yml will provide better restrictions on installation.

edge problem? I'm getting vertical of zero pixels at tile edge here

I want to use the glo30 DEM with ISCE. I basically use the notebook here, 'staging for isce2', but with a boundary [-74, 41,-70, 44].
The result is a vertical line of 0 valued pixels (~1 pixel wide) going across the stitched tile at about the -72.00000 coordinate through the entire region. This results in a problem with ISCE2 processing. Here is a picture:

Software Versioning, ASF CI/CD

• automate software versioning with github actions (as here including labeled-pr, changelog, and test-and-tag)
• Automate pypi release as here

This is now done with s1-enumerator and will be replicated.

Removal of `2021` from AWS Glo Registry

The bug

The example on the github homepage does not work due to removal of 2021 prefix in glo-30 aws bucket.

To Reproduce

bounds = [-119.085, 33.402, -118.984, 35.435]
X, p = stitch_dem(bounds,
                  dem_name='glo_30',
                  dst_ellipsoidal_height=False,
                  dst_area_or_point='Point')

---------------------------------------------------------------------------
CPLE_HttpResponseError                    Traceback (most recent call last)
File rasterio/_base.pyx:261, in rasterio._base.DatasetBase.__init__()

File rasterio/_shim.pyx:78, in rasterio._shim.open_dataset()

File rasterio/_err.pyx:216, in rasterio._err.exc_wrap_pointer()

CPLE_HttpResponseError: HTTP response code: 404

During handling of the above exception, another exception occurred:

RasterioIOError                           Traceback (most recent call last)
Input In [22], in <cell line: 2>()
      1 bounds = [-119.085, 33.402, -118.984, 35.435]
----> 2 X, p = stitch_dem(bounds,
      3                   dem_name='glo_30',
      4                   dst_ellipsoidal_height=False,
      5                   dst_area_or_point='Point')

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/site-packages/dem_stitcher/stitcher.py:304, in stitch_dem(bounds, dem_name, dst_ellipsoidal_height, dst_area_or_point, dst_resolution, n_threads_reproj, n_threads_downloading, driver, fill_in_glo_30)
    302     return RASTER_READERS[dem_name](url)
    303 with ThreadPoolExecutor(max_workers=n_threads_downloading) as executor:
--> 304     results = list(tqdm(executor.map(reader, urls),
    305                         total=len(urls),
    306                         desc=f'Reading {dem_name} Datasets'))
    308 # If datasets are non-existent, returns None
    309 datasets = list(filter(lambda x: x is not None, results))

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/site-packages/tqdm/std.py:1195, in tqdm.__iter__(self)
   1192 time = self._time
   1194 try:
-> 1195     for obj in iterable:
   1196         yield obj
   1197         # Update and possibly print the progressbar.
   1198         # Note: does not call self.update(1) for speed optimisation.

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/concurrent/futures/_base.py:609, in Executor.map.<locals>.result_iterator()
    606 while fs:
    607     # Careful not to keep a reference to the popped future
    608     if timeout is None:
--> 609         yield fs.pop().result()
    610     else:
    611         yield fs.pop().result(end_time - time.monotonic())

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/concurrent/futures/_base.py:446, in Future.result(self, timeout)
    444     raise CancelledError()
    445 elif self._state == FINISHED:
--> 446     return self.__get_result()
    447 else:
    448     raise TimeoutError()

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/concurrent/futures/_base.py:391, in Future.__get_result(self)
    389 if self._exception:
    390     try:
--> 391         raise self._exception
    392     finally:
    393         # Break a reference cycle with the exception in self._exception
    394         self = None

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/concurrent/futures/thread.py:58, in _WorkItem.run(self)
     55     return
     57 try:
---> 58     result = self.fn(*self.args, **self.kwargs)
     59 except BaseException as exc:
     60     self.future.set_exception(exc)

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/site-packages/dem_stitcher/stitcher.py:302, in stitch_dem.<locals>.reader(url)
    301 def reader(url):
--> 302     return RASTER_READERS[dem_name](url)

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/site-packages/dem_stitcher/dem_readers.py:12, in read_dem(dem_path)
     11 def read_dem(dem_path: str) -> rasterio.DatasetReader:
---> 12     ds = rasterio.open(dem_path)
     13     return ds

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/site-packages/rasterio/env.py:437, in ensure_env_with_credentials.<locals>.wrapper(*args, **kwds)
    434     session = DummySession()
    436 with env_ctor(session=session):
--> 437     return f(*args, **kwds)

File ~/opt/anaconda3/envs/dem-stitcher/lib/python3.9/site-packages/rasterio/__init__.py:220, in open(fp, mode, driver, width, height, count, crs, transform, dtype, nodata, sharing, **kwargs)
    216 # Create dataset instances and pass the given env, which will
    217 # be taken over by the dataset's context manager if it is not
    218 # None.
    219 if mode == 'r':
--> 220     s = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
    221 elif mode == "r+":
    222     s = get_writer_for_path(path, driver=driver)(
    223         path, mode, driver=driver, sharing=sharing, **kwargs
    224     )

File rasterio/_base.pyx:263, in rasterio._base.DatasetBase.__init__()

RasterioIOError: HTTP response code: 404

Missing glo-30 DEM tiles over Azerbaijan, Armenia

The bug

Glo-30 tiles have missing tiles over said areas. There are downstream implications related to removing topographic phase for GUNWs.

A reference to some of the missing tiles: https://twitter.com/EricFielding/status/1531639815610830848

Country identifier | GeoCell ID
Armenia | N38E046
Armenia | N39E044
Armenia | N39E045
Armenia | N39E046
Armenia | N40E043
Armenia | N40E044
Armenia | N40E045
Armenia | N41E043
Armenia | N41E044
Armenia | N41E045
Azerbaijan | N38E046
Azerbaijan | N39E045
Azerbaijan | N39E046
Azerbaijan | N39E047
Azerbaijan | N40E045
Azerbaijan | N40E046
Azerbaijan | N40E047
Azerbaijan | N38E045
Azerbaijan | N38E046
Azerbaijan | N38E048
Azerbaijan | N38E049
Azerbaijan | N39E044
Azerbaijan | N39E045
Azerbaijan | N39E046
Azerbaijan | N39E047
Azerbaijan | N39E048
Azerbaijan | N39E049
Azerbaijan | N40E044
Azerbaijan | N40E045
Azerbaijan | N40E046
Azerbaijan | N40E047
Azerbaijan | N40E048
Azerbaijan | N40E049
Azerbaijan | N40E050
Azerbaijan | N41E044
Azerbaijan | N41E045
Azerbaijan | N41E046
Azerbaijan | N41E047
Azerbaijan | N41E048
Azerbaijan | N41E049

To Reproduce

Check this package data file and verify above tiles doe not exist

Additional context

Downstream topographic phase removal - see Issue #67

Access NASADEM Water Body layer (SWB)

Hi!

First of all, thanks for the great package. I was wondering if there is a simple way to access the water bodies layer of the NASADEM-HGT product. It is listed here as the third layer. I assume your stitch_dem function always returns the first layer only?

In the merge_and_transform_dem_tiles function there are two calls where you set:

dem_arr = dem_arr[0, ...]

I assume this is where you drop all but the first layer (i.e. all but the DEM layer)?

I hope the question is somewhat clear; please let me know if not.

Best,
Konstantin

Verification of DEM heights and Geoid removal in comparison to independent dataset

Ensure correctness of geoid removal using ESA Copernicus data as indicated here.

Anti-meridian

The bug

At the antimeridian line, there is no "unwrapping" of the tiles and therefore only longitudes -180 to 180 are included.

To Reproduce

bounds = [-180.25, 51.25, -179.75, 51.75]
X, p = stitch_dem(bounds,
                               dem_name='glo_30')

Additional context

The above should yield the same as

bounds = [179.75, 51.25, 180.25, 51.75]
X, p = stitch_dem(bounds,
                               dem_name='glo_30')

Update ISCE2 Ingest and Clarification of how to stage DEMs

Do not want to include ISCE2 as requirement here.

Will clarify how to stage DEMs for ISCE2 including using fixImageXML.py from ISCE2 as done here.

Basically, make sure:

1. driver is ISCE
2. Tagged with ellipsoidal heights
3. Run the paths in the XML metadata file.

Left, bottom boundaries of the DEM are getting rounded to an integer

The bug

While the right/top boundaries are exactly as specified, the left/bottm (western/southern) boundaries seem to get rounded.

To Reproduce

(following the README)

import rasterio
bounds = [176.678207, 50.908962, 179.697601, 52.754662]
X, p = dem_stitcher.stitch_dem(bounds, dem_name='glo_30')
with rasterio.open('dem.tif', 'w', **p) as ds:
    ds.write(X2, 1)
    ds.update_tags(AREA_OR_POINT='Point')

$ rio bounds --bbox dem.tif
[177.0, 51.0, 179.69791666666666, 52.75472222222222]

Additional context

In [36]: rio.show_versions()
rasterio info:
  rasterio: 1.3.8
      GDAL: 3.7.1
      PROJ: 9.2.1
      GEOS: 3.12.0
 PROJ DATA: /Users/staniewi/miniconda3/envs/dem-env/share/proj
 GDAL DATA: /Users/staniewi/miniconda3/envs/dem-env/share/gdal

System:
    python: 3.11.5 | packaged by conda-forge | (main, Aug 27 2023, 03:33:12) [Clang 15.0.7 ]
executable: /Users/staniewi/miniconda3/envs/dem-env/bin/python3.11
   machine: macOS-13.5.1-arm64-arm-64bit

Python deps:
    affine: 2.4.0
     attrs: 23.1.0
   certifi: 2023.07.22
     click: 8.1.7
     cligj: 0.7.2
    cython: None
     numpy: 1.25.2
    snuggs: 1.4.7
click-plugins: None
setuptools: 68.1.2