Remove pygmt library usage because the recent version crashes Jupyter…

… kernel in Docker containers
AlexeyPechnikov · Feb 26, 2024 · 7281edf · 7281edf
1 parent de11de6
commit 7281edf
Show file tree

Hide file tree

Showing 5 changed files with 181 additions and 122 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -60,7 +60,6 @@ def get_version():
                       'h5py',
                       'nc-time-axis',
                       'statsmodels>=0.14.0',
-                      'pygmt',
                       'remotezip',
                       'asf_search',
                       'imageio',