From 7281edf90ad471c03c6e2ff0a87beb361e48cc4e Mon Sep 17 00:00:00 2001 From: Alexey Pechnikov Date: Mon, 26 Feb 2024 16:03:52 +0700 Subject: [PATCH] Remove pygmt library usage because the recent version crashes Jupyter kernel in Docker containers --- pygmtsar/pygmtsar/GMT.py | 218 ++++++++++++++++++++++------ pygmtsar/pygmtsar/Stack_dem.py | 17 +-- pygmtsar/pygmtsar/Stack_landmask.py | 65 +-------- pygmtsar/pygmtsar/__init__.py | 2 +- pygmtsar/setup.py | 1 - 5 files changed, 181 insertions(+), 122 deletions(-) diff --git a/pygmtsar/pygmtsar/GMT.py b/pygmtsar/pygmtsar/GMT.py index 46a6cbef..c8302d20 100644 --- a/pygmtsar/pygmtsar/GMT.py +++ b/pygmtsar/pygmtsar/GMT.py @@ -11,9 +11,121 @@ from .tqdm_joblib import tqdm_joblib class GMT(datagrid, tqdm_joblib): - - # https://docs.generic-mapping-tools.org/6.0/datasets/earth_relief.html - # only bicubic interpolation supported as the best one for the case + +# 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 + def download_dem(self, geometry, filename=None, product='1s', skip_exist=True): """ Download and preprocess SRTM digital elevation model (DEM) data using GMT library. @@ -35,52 +147,62 @@ def download_dem(self, geometry, filename=None, product='1s', skip_exist=True): 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') + + Notes + -------- + https://docs.generic-mapping-tools.org/6.0/datasets/earth_relief.html """ - import numpy as np - import pygmt - # suppress warnings - pygmt.config(GMT_VERBOSE='errors') 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: + print(f'Error executing gmt grdcut: {stderr_data}') + return stdout_data.strip() 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 + # 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 - 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]) + 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 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 + if filename is None: + da = xr.open_dataarray(grdname, engine=self.netcdf_engine, chunks=self.netcdf_chunksize) + os.remove(grdname) + return da - def download_landmask(self, geometry, filename=None, product='1s', skip_exist=True): + def download_landmask(self, geometry, filename=None, product='1s', resolution='f', skip_exist=True): """ Download the landmask and save as NetCDF file. @@ -98,14 +220,23 @@ def download_landmask(self, geometry, filename=None, product='1s', skip_exist=Tr ----- 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 + 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}', '-N1/NaN', 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: + print(f'Error executing grdlandmask: {stderr_data}') + return stdout_data.strip() 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') @@ -115,15 +246,16 @@ def download_landmask(self, geometry, filename=None, product='1s', skip_exist=Tr 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) + 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: - landmask = pygmt.grdlandmask(resolution='f', region=[lon_start, lon_end, lat_start, lat_end], spacing=product, maskvalues='NaN/1') + grdlandmask(self.get_bounds(geometry), product, resolution, grdname) 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 + if filename is None: + da = xr.open_dataarray(grdname, engine=self.netcdf_engine, chunks=self.netcdf_chunksize) + os.remove(grdname) + return da diff --git a/pygmtsar/pygmtsar/Stack_dem.py b/pygmtsar/pygmtsar/Stack_dem.py index 018a5d7c..3bad02b7 100644 --- a/pygmtsar/pygmtsar/Stack_dem.py +++ b/pygmtsar/pygmtsar/Stack_dem.py @@ -27,11 +27,9 @@ def get_extent_ra(self): def get_extent(self, grid=None, subswath=None): import numpy as np - extent = self.get_reference(subswath).dissolve().envelope.item() + bounds = self.get_bounds(self.get_reference(subswath)) if grid is None: - return extent - bounds = np.round(extent.bounds, 3) - #print ('xmin, xmax', xmin, xmax) + return bounds return grid\ .transpose('lat','lon')\ .sel(lat=slice(bounds[1], bounds[3]), @@ -226,17 +224,6 @@ def load_dem(self, data, geometry='auto'): print ('ERROR: argument is not an Xarray object and it is not a file name') # crop -# if type(geometry) == str and geometry == 'auto': -# # apply scenes geometry -# extent = self.get_extent() -# elif isinstance(geometry, gpd.GeoDataFrame): -# extent = geometry.dissolve().envelope.item() -# elif isinstance(geometry, gpd.GeoSeries): -# geometry = geometry.unary_union.envelope -# # round the coordinates up to 1m -# #minx, miny, maxx, maxy = np.round(geometry.bounds, 5) -# #print ('minx, miny, maxx, maxy', minx, miny, maxx, maxy) -# bounds = np.round(extent.bounds, 5) bounds = self.get_bounds(self.get_extent() if type(geometry) == str and geometry == 'auto' else geometry) ortho = ortho\ .transpose('lat','lon')\ diff --git a/pygmtsar/pygmtsar/Stack_landmask.py b/pygmtsar/pygmtsar/Stack_landmask.py index ce1cc403..b0d9218f 100644 --- a/pygmtsar/pygmtsar/Stack_landmask.py +++ b/pygmtsar/pygmtsar/Stack_landmask.py @@ -81,56 +81,8 @@ def get_landmask(self): return landmask def download_landmask(self, product='1s', debug=False): - print ('NOTE: Function is deprecated. Download land mask using GMT.download_landmask(stack.get_dem())') - print ('function and load with stack.load_landmask() function.') - - """ - 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). - debug : bool, optional - Whether to enable debug mode. Default is False. - - Examples - -------- - stack.download_landmask() - - Notes - ----- - This method downloads the landmask using GMT's local data or server. The landmask is built based on the interferogram DEM area. - """ - import pygmt - import os - from tqdm.auto import tqdm - - if self.landmask_filename is not None: - print ('NOTE: landmask exists, ignore the command. Use Stack.set_landmask(None) to allow new landmask downloading') - return - - # generate the same as DEM grid - landmask_filename = os.path.join(self.basedir, 'landmask.nc') - - # geographical grid for interferogram area only - dem = self.get_dem() - llmin = dem.lon.min().item() - llmax = dem.lon.max().item() - ltmin = dem.lat.min().item() - ltmax = dem.lat.max().item() - region = f'{llmin}/{llmax}/{ltmin}/{ltmax}' - if debug: - print('region', region) - - with tqdm(desc='Landmask Downloading', total=1) as pbar: - landmask = pygmt.grdlandmask(resolution='f', region=region, spacing=product, maskvalues='NaN/1') - if os.path.exists(landmask_filename): - os.remove(landmask_filename) - landmask.to_netcdf(landmask_filename) - pbar.update(1) - - self.landmask_filename = landmask_filename + print ('NOTE: Function is deprecated. Download land mask using GMT.download_landmask()') + print ('function and load with Stack.load_landmask() function.') def load_landmask(self, data, geometry='auto'): """ @@ -188,18 +140,7 @@ def load_landmask(self, data, geometry='auto'): else: print ('ERROR: argument is not an Xarray object and it is not a file name') -# # crop -# if type(geometry) == str and geometry == 'auto': -# # apply scenes geometry -# extent = self.get_extent() -# elif isinstance(geometry, gpd.GeoDataFrame): -# extent = geometry.dissolve().envelope.item() -# elif isinstance(geometry, gpd.GeoSeries): -# geometry = geometry.unary_union.envelope -# # round the coordinates up to 1m -# #minx, miny, maxx, maxy = np.round(geometry.bounds, 5) -# #print ('minx, miny, maxx, maxy', minx, miny, maxx, maxy) -# bounds = np.round(extent.bounds, 5) + # crop bounds = self.get_bounds(self.get_extent() if type(geometry) == str and geometry == 'auto' else geometry) landmask = landmask\ .transpose('lat','lon')\ diff --git a/pygmtsar/pygmtsar/__init__.py b/pygmtsar/pygmtsar/__init__.py index b20025ce..68d3808d 100644 --- a/pygmtsar/pygmtsar/__init__.py +++ b/pygmtsar/pygmtsar/__init__.py @@ -7,7 +7,7 @@ # # Licensed under the BSD 3-Clause License (see LICENSE for details) # ---------------------------------------------------------------------------- -__version__ = '2024.2.21.post1' +__version__ = '2024.2.21.post3' # unified progress indicators from .tqdm_joblib import tqdm_joblib diff --git a/pygmtsar/setup.py b/pygmtsar/setup.py index 72185d2f..1223a055 100644 --- a/pygmtsar/setup.py +++ b/pygmtsar/setup.py @@ -60,7 +60,6 @@ def get_version(): 'h5py', 'nc-time-axis', 'statsmodels>=0.14.0', - 'pygmt', 'remotezip', 'asf_search', 'imageio',