ACCESS-Cloud-Based-InSAR / dem-stitcher

Download and merge DEM tiles
Apache License 2.0
39 stars 15 forks source link

Access NASADEM Water Body layer (SWB) #55

Closed konstantinklemmer closed 1 year ago

konstantinklemmer commented 1 year ago

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

cmarshak commented 1 year ago

Hey Konstantin,

Glad it's helpful.

Unfortunately, it's not in the scope of this package as is. I can share with you the rough idea of how you might approach this.

  1. The urls will be the same as for NASADEM.
from dem_stitcher.datasets import get_overlapping_dem_tiles

urls = get_overlapping_dem_tiles([-121.5, 34.95, -120.2, 36.25], 'nasadem').url.tolist()

Note: if you unzip a url, you will find the *.swb layer in the unzipped directory. Here is an example url: https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_HGT.001/2000.02.11/NASADEM_HGT_n34w121.zip

  1. Modify these mostly stand-alone functions (or reorganize them) to read the swb layer instead of the hgt layer. Might need a little bit of work (haven't tried it) because SRTM has its own gdal driver: https://gdal.org/drivers/raster/srtmhgt.html, but swb does not. You could infer the geometadata from the hgt file and then just generate a dataset.

Hope this helps. Good luck!

konstantinklemmer commented 1 year ago

Thanks for the swift reply!

cmarshak commented 1 year ago

Here is another option from @mgovorcin using ISCE2:

import isce
from isceobj.Alos2Proc.runDownloadDem import download_wbd
# Download SRTM-SWDB water mask
# ISCE2 requires South, North, East, West extents (or ymin, ymax, xmin, xmax)
mask_filename = download_wbd(extent[1], extent[3], extent[0], extent[2])