From 94815de8d0338b84c240d2a9386f8e28ed93d764 Mon Sep 17 00:00:00 2001 From: Alexey Pechnikov Date: Thu, 14 Mar 2024 13:51:01 +0700 Subject: [PATCH] Prepare to remove modules AWS and GMT to use Tiles module instead --- pygmtsar/pygmtsar/AWS.py | 168 +------------------------------ pygmtsar/pygmtsar/GMT.py | 208 +-------------------------------------- 2 files changed, 5 insertions(+), 371 deletions(-) diff --git a/pygmtsar/pygmtsar/AWS.py b/pygmtsar/pygmtsar/AWS.py index 6ed4e94d..8d4cfd29 100644 --- a/pygmtsar/pygmtsar/AWS.py +++ b/pygmtsar/pygmtsar/AWS.py @@ -7,126 +7,7 @@ # # Licensed under the BSD 3-Clause License (see LICENSE for details) # ---------------------------------------------------------------------------- -from .datagrid import datagrid -from .tqdm_joblib import tqdm_joblib - -class AWS(datagrid, tqdm_joblib): - # https://copernicus-dem-90m.s3.eu-central-1.amazonaws.com - # https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N38_00_E038_00_DEM/Copernicus_DSM_COG_10_N38_00_E038_00_DEM.tif - base_url_glo = 'https://copernicus-dem-{resolution}m.s3.amazonaws.com' - path_id_glo = 'Copernicus_DSM_COG_{product1}0_{SN2}_00_{WE3}_00_DEM' - # Copernicus_DSM_COG_10_S16_00_W043_00_DEM - tile_id_glo = path_id_glo + '.tif' - - #aws s3 ls --no-sign-request s3://elevation-tiles-prod/skadi/ - # https://s3.amazonaws.com/elevation-tiles-prod/skadi/N20/N20E000.hgt.gz - base_url_srtm = 'https://s3.amazonaws.com/elevation-tiles-prod/skadi' - path_id_srtm = '{SN2}' - tile_id_srtm = '{SN2}{WE3}.hgt.gz' - - def _download_tile_glo(self, product, lon, lat, debug=False): - """ - Download Copernicus GLO-30/GLO-90 Digital Elevation Model tiles from open AWS storage. - - Previously, it was in-memory processing but it seems not stable on slow internet connection: - import io - with io.BytesIO(response.content) as f: - #tile = xr.open_dataarray(f, engine='rasterio')... - tile = rio.open_rasterio(f, chunks=self.chunksize)... - """ - import rioxarray as rio - import requests - import os - import tempfile - - product1 = int(product[0]) - SN2 = f'{"S" if lat<0 else "N"}{abs(lat):02}' - WE3 = f'{"W" if lon<0 else "E"}{abs(lon):03}' - base_url = self.base_url_glo.format(resolution=30 if product=='1s' else 90) - path_id = self.path_id_glo.format(product1=product1, SN2=SN2, WE3=WE3) - tile_id = self.tile_id_glo.format(product1=product1, SN2=SN2, WE3=WE3) - tile_url = f'{base_url}/{path_id}/{tile_id}' - tile_filename = os.path.join(tempfile.gettempdir(), tile_id) - if debug: - print ('DEBUG _download_tile_glo: tile_url', tile_url) - print ('DEBUG _download_tile_glo: tile_filename', tile_filename) - try: - with requests.get(tile_url, stream=True) as response: - response.raise_for_status() - with open(tile_filename, 'wb') as f: - for chunk in response.iter_content(chunk_size=256*1024): - f.write(chunk) - with rio.open_rasterio(tile_filename) as raster: - tile = raster\ - .squeeze(drop=True)\ - .rename({'y': 'lat', 'x': 'lon'})\ - .drop_vars('spatial_ref')\ - .load() - if tile.lat.diff('lat')[0].item() < 0: - tile = tile.reindex(lat=tile.lat[::-1]) - except requests.exceptions.RequestException as e: - # offshore tiles are missed by design - print(f'Request error for {tile_id}: {e}') - tile = None - except Exception as e: - print(e) - raise - finally: - if os.path.exists(tile_filename): - os.remove(tile_filename) - return tile - - def _download_tile_srtm(self, product, lon, lat, debug=False): - """ - Download NASA SRTM Digital Elevation Model tiles from open AWS storage. - """ - import rioxarray as rio - import requests - import os - import tempfile - import gzip - - SN2 = f'{"S" if lat<0 else "N"}{abs(lat):02}' - WE3 = f'{"W" if lon<0 else "E"}{abs(lon):03}' - base_url = self.base_url_srtm - path_id = self.path_id_srtm.format(SN2=SN2, WE3=WE3) - tile_id = self.tile_id_srtm.format(SN2=SN2, WE3=WE3) - tile_url = f'{base_url}/{path_id}/{tile_id}' - # remove .gz extension - tile_filename = os.path.join(tempfile.gettempdir(), tile_id[:-3]) - if debug: - print ('DEBUG _download_tile_srtm: tile_url', tile_url) - print ('DEBUG _download_tile_srtm: tile_filename', tile_filename) - try: - with requests.get(tile_url, stream=True) as response: - response.raise_for_status() - # Stream the content under the context of gzip decompression - with gzip.GzipFile(fileobj=response.raw) as gz: - with open(tile_filename, 'wb') as f: - while True: - chunk = gz.read(256 * 1024) - if not chunk: - break - f.write(chunk) - with rio.open_rasterio(tile_filename) as raster: - tile = raster\ - .squeeze(drop=True)\ - .rename({'y': 'lat', 'x': 'lon'})\ - .drop_vars('spatial_ref')\ - .load() - if tile.lat.diff('lat')[0].item() < 0: - tile = tile.reindex(lat=tile.lat[::-1]) - except requests.exceptions.RequestException as e: - # offshore tiles are missed by design - print(f'Request error for {tile_id}: {e}') - tile = None - except Exception as e: - print(e) - raise - finally: - if os.path.exists(tile_filename): - os.remove(tile_filename) - return tile +class AWS(): def download_dem(self, geometry, filename=None, product='1s', provider='GLO', n_jobs=4, joblib_backend='loky', skip_exist=True): """ @@ -138,49 +19,4 @@ def download_dem(self, geometry, filename=None, product='1s', provider='GLO', n_ AWS().download_dem(S1.scan_slc(DATADIR), 'dem.nc') """ - import xarray as xr - import numpy as np - from tqdm.auto import tqdm - import joblib - import os - - assert provider in ['GLO', 'SRTM'], f'ERROR: provider name is invalid: {provider}. Expected names are "GLO", "SRTM".' - if provider == 'SRTM': - assert product in ['1s'], f'ERROR: only product="1s" is supported for provider {provider}.' - if provider == 'GLO': - assert product in ['1s', '3s'], f'ERROR: product name is invalid: {product} for provider {provider}. Expected names are "1s", "3s".' - - if filename is not None and os.path.exists(filename) and skip_exist: - print ('NOTE: DEM file exists, ignore the command. Use "skip_exist=False" or omit the filename to allow new downloading.') - return - - bounds = self.get_bounds(geometry) - - # it produces 4 tiles for cases like (39.5, 39.5, 40.0, 40.0) - #left, right = int(np.floor(lon_start)), int(np.floor(lon_end)) - #bottom, top = int(np.floor(lat_start)), int(np.floor(lat_end)) - # enhancement to produce a single tile for cases like (39.5, 39.5, 40.0, 40.0) - lon_start, lat_start, lon_end, lat_end = bounds - left = np.floor(min(lon_start, lon_end)) - right = np.ceil(max(lon_start, lon_end)) - 1 - bottom = np.floor(min(lat_start, lat_end)) - top = np.ceil(max(lat_start, lat_end)) - 1 - left, right = int(left), int(right) - bottom, top = int(bottom), int(top) - #print ('left, right', left, right, 'bottom, top', bottom, top) - - job_tile = self._download_tile_glo if provider == 'GLO' else self._download_tile_srtm - with self.tqdm_joblib(tqdm(desc=f'AWS {provider} DEM Tiles Downloading', total=(right-left+1)*(top-bottom+1))) as progress_bar: - tile_xarrays = joblib.Parallel(n_jobs=n_jobs, backend=joblib_backend)(joblib.delayed(job_tile)(product, x, y)\ - for x in range(left, right + 1) for y in range(bottom, top + 1)) - - dem = xr.combine_by_coords([tile for tile in tile_xarrays if tile is not None]) - dem = dem.sel(lat=slice(bounds[1], bounds[3]), lon=slice(bounds[0], bounds[2])) - - if filename is not None: - if os.path.exists(filename): - os.remove(filename) - encoding = {'dem': self._compression(dem.shape)} - dem.rename('dem').to_netcdf(filename, encoding=encoding, engine=self.netcdf_engine) - else: - return dem + print('NOTE: AWS module is removed. Download DEM using Tiles().download_dem(AOI, filename=...).') diff --git a/pygmtsar/pygmtsar/GMT.py b/pygmtsar/pygmtsar/GMT.py index 24f5839a..b118ac77 100644 --- a/pygmtsar/pygmtsar/GMT.py +++ b/pygmtsar/pygmtsar/GMT.py @@ -7,125 +7,7 @@ # # Licensed under the BSD 3-Clause License (see LICENSE for details) # ---------------------------------------------------------------------------- -from .datagrid import datagrid -from .tqdm_joblib import tqdm_joblib - -class GMT(datagrid, tqdm_joblib): - -# def download_dem(self, geometry, filename=None, product='1s', skip_exist=True): -# """ -# Download and preprocess SRTM digital elevation model (DEM) data using GMT library. -# -# Parameters -# ---------- -# product : str, optional -# Product type of the DEM data. Available options are '1s' or 'SRTM1' (1 arcsec ~= 30m, default) -# and '3s' or 'SRTM3' (3 arcsec ~= 90m). -# -# Returns -# ------- -# None or Xarray Dataarray -# -# Examples -# -------- -# Download default STRM1 DEM (~30 meters): -# GMT().download_dem() -# -# Download STRM3 DEM (~90 meters): -# GMT.download_dem(product='STRM3') -# -# Download default STRM DEM to cover the selected area AOI: -# GMT().download_dem(AOI) -# -# Download default STRM DEM to cover all the scenes: -# GMT().download_dem(stack.get_extent().buffer(0.1)) -# -# import pygmt -# # Set the GMT data server limit to N Mb to allow for remote downloads -# pygmt.config(GMT_DATA_SERVER_LIMIT=1e6) -# GMT().download_dem(AOI, product='1s') -# """ -# import numpy as np -# import pygmt -# # suppress warnings -# pygmt.config(GMT_VERBOSE='errors') -# import os -# from tqdm.auto import tqdm -# -# if product in ['SRTM1', '1s', '01s']: -# resolution = '01s' -# elif product in ['SRTM3', '3s', '03s']: -# resolution = '03s' -# else: -# print (f'ERROR: unknown product {product}. Available only SRTM1 ("01s") and SRTM3 ("03s") DEM using GMT servers') -# return -# -# if filename is not None and os.path.exists(filename) and skip_exist: -# print ('NOTE: DEM file exists, ignore the command. Use "skip_exist=False" or omit the filename to allow new downloading') -# return -# -# lon_start, lat_start, lon_end, lat_end = self.get_bounds(geometry) -# with tqdm(desc='GMT SRTM DEM Downloading', total=1) as pbar: -# # download DEM using GMT extent W E S N -# dem = pygmt.datasets.load_earth_relief(resolution=resolution, region=[lon_start, lon_end, lat_start, lat_end]) -# pbar.update(1) -# -# if filename is not None: -# if os.path.exists(filename): -# os.remove(filename) -# encoding = {'dem': self._compression(dem.shape)} -# dem.rename('dem').load().to_netcdf(filename, encoding=encoding, engine=self.netcdf_engine) -# else: -# return dem - -# def download_landmask(self, geometry, filename=None, product='1s', resolution='f', skip_exist=True): -# """ -# Download the landmask and save as NetCDF file. -# -# Parameters -# ---------- -# product : str, optional -# Available options are '1s' (1 arcsec ~= 30m, default) and '3s' (3 arcsec ~= 90m). -# -# Examples -# -------- -# from pygmtsar import GMT -# landmask = GMT().download_landmask(stack.get_dem()) -# -# Notes -# ----- -# This method downloads the landmask using GMT's local data or server. -# """ -# import geopandas as gpd -# import numpy as np -# import pygmt -# # suppress warnings -# pygmt.config(GMT_VERBOSE='errors') -# import os -# #import subprocess -# from tqdm.auto import tqdm -# -# if filename is not None and os.path.exists(filename) and skip_exist: -# print ('NOTE: landmask file exists, ignore the command. Use "skip_exist=False" or omit the filename to allow new downloading') -# return -# -# if not product in ['1s', '3s']: -# print (f'ERROR: unknown product {product}. Available only "1s" or "3s" land mask using GMT servers') -# return -# -# lon_start, lat_start, lon_end, lat_end = self.get_bounds(geometry) -# with tqdm(desc='GMT Landmask Downloading', total=1) as pbar: -# landmask = pygmt.grdlandmask(resolution=resolution, region=[lon_start, lon_end, lat_start, lat_end], spacing=product, maskvalues='NaN/1') -# pbar.update(1) -# -# if filename is not None: -# if os.path.exists(filename): -# os.remove(filename) -# encoding = {'landmask': self._compression(landmask.shape)} -# landmask.rename('landmask').load().to_netcdf(filename, encoding=encoding, engine=self.netcdf_engine) -# else: -# return landmask - +class GMT(): def download_dem(self, geometry, filename=None, product='1s', skip_exist=True): """ Download and preprocess SRTM digital elevation model (DEM) data using GMT library. @@ -158,51 +40,7 @@ def download_dem(self, geometry, filename=None, product='1s', skip_exist=True): -------- https://docs.generic-mapping-tools.org/6.0/datasets/earth_relief.html """ - import xarray as xr - import os - import subprocess - from tqdm.auto import tqdm - import tempfile - import warnings - warnings.filterwarnings('ignore') - - def load_earth_relief(bounds, product, filename): - if os.path.exists(filename): - os.remove(filename) - argv = ['gmt', 'grdcut', f'@earth_relief_{product}', f'-R{bounds[0]}/{bounds[2]}/{bounds[1]}/{bounds[3]}', f'-G{filename}'] - #print ('gmt grdcut argv:', ' '.join(argv)) - p = subprocess.Popen(argv, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf8') - stdout_data, stderr_data = p.communicate() - if p.returncode != 0 or not os.path.exists(filename): - print('Error executing external command "gmt grdcut":') - raise ValueError(stderr_data) - return stdout_data - - if product in ['SRTM1', '1s', '01s']: - resolution = '01s' - elif product in ['SRTM3', '3s', '03s']: - resolution = '03s' - else: - # use the specified string as is - pass - - if filename is not None and os.path.exists(filename) and skip_exist: - print ('NOTE: DEM file exists, ignore the command. Use "skip_exist=False" or omit the filename to allow new downloading') - return - - grdname = filename - if grdname is None: - with tempfile.NamedTemporaryFile() as tmpfile: - grdname = tmpfile.name + '.grd' - - with tqdm(desc=f'GMT SRTM {product} DEM Downloading', total=1) as pbar: - load_earth_relief(self.get_bounds(geometry), resolution, grdname) - pbar.update(1) - - if filename is None: - da = xr.open_dataarray(grdname, engine=self.netcdf_engine).load() - os.remove(grdname) - return da + print('NOTE: GMT module is removed. Download DEM using Tiles().download_dem(AOI, filename=...).') def download_landmask(self, geometry, filename=None, product='1s', resolution='f', skip_exist=True): """ @@ -222,44 +60,4 @@ def download_landmask(self, geometry, filename=None, product='1s', resolution='f ----- This method downloads the landmask using GMT's local data or server. """ - import xarray as xr - import os - import subprocess - from tqdm.auto import tqdm - import tempfile - import warnings - warnings.filterwarnings('ignore') - - def grdlandmask(bounds, product, resolution, filename): - if os.path.exists(filename): - os.remove(filename) - argv = ['gmt', 'grdlandmask', f'-R{bounds[0]}/{bounds[2]}/{bounds[1]}/{bounds[3]}', f'-I{product}', f'-D{resolution}', '-NNaN/1', '-rp', f'-G{filename}'] - #print ('grdlandmask argv:', ' '.join(argv)) - p = subprocess.Popen(argv, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf8') - stdout_data, stderr_data = p.communicate() - if p.returncode != 0 or not os.path.exists(filename): - print('Error executing external command "gmt grdlandmask":') - raise ValueError(stderr_data) - return stdout_data - - if filename is not None and os.path.exists(filename) and skip_exist: - print ('NOTE: landmask file exists, ignore the command. Use "skip_exist=False" or omit the filename to allow new downloading') - return - - if not product in ['1s', '3s']: - print (f'ERROR: unknown product {product}. Available only "1s" or "3s" land mask using GMT servers') - return - - grdname = filename - if grdname is None: - with tempfile.NamedTemporaryFile() as tmpfile: - grdname = tmpfile.name + '.grd' - - with tqdm(desc='GMT Landmask Downloading', total=1) as pbar: - grdlandmask(self.get_bounds(geometry), product, resolution, grdname) - pbar.update(1) - - if filename is None: - da = xr.open_dataarray(grdname, engine=self.netcdf_engine).load() - os.remove(grdname) - return da + print('NOTE: GMT module is removed. Download landmask using Tiles().download_landmask(AOI, filename=...).')