From 50bbe59f87a4b3750209395e80f4b40bbd114f53 Mon Sep 17 00:00:00 2001 From: Alexey Pechnikov Date: Tue, 8 Aug 2023 20:30:50 +0700 Subject: [PATCH] Allow to use local DEM file for function download_dem() --- pygmtsar/pygmtsar/SBAS_dem.py | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/pygmtsar/pygmtsar/SBAS_dem.py b/pygmtsar/pygmtsar/SBAS_dem.py index efd2468c..a65db705 100644 --- a/pygmtsar/pygmtsar/SBAS_dem.py +++ b/pygmtsar/pygmtsar/SBAS_dem.py @@ -121,6 +121,7 @@ def get_dem(self, subswath=None, geoloc=False, buffer_degrees=0.02): .sel(lat=slice(bounds[1]-buffer_degrees, bounds[3]+buffer_degrees), lon=slice(bounds[0]-buffer_degrees, bounds[2]+buffer_degrees)) + # buffer required to get correct (binary) results from SAT_llt2rat tool # small margin produces insufficient DEM not covers the defined area # https://docs.generic-mapping-tools.org/6.0/datasets/earth_relief.html @@ -159,6 +160,12 @@ def download_dem(self, backend=None, product='SRTM1', resolution_meters=None, me Download STRM3 DEM with a resolution of 90 meters and convert it to a 120-meter grid: sbas.download_dem(product='STRM3', resolution_meters=120) + Load and crop from local NetCDF file: + sbas.download_dem(product='GEBCO_2020/GEBCO_2020.nc') + + Load and crop from local GeoTIF file: + sbas.download_dem(product='GEBCO_2019.tif') + Notes ----- This method uses the GMT servers to download SRTM 1 or 3 arc-second DEM data. The downloaded data is then preprocessed @@ -168,6 +175,7 @@ def download_dem(self, backend=None, product='SRTM1', resolution_meters=None, me import xarray as xr import numpy as np import pygmt + import rioxarray as rio import os #import subprocess from tqdm.auto import tqdm @@ -190,7 +198,9 @@ def download_dem(self, backend=None, product='SRTM1', resolution_meters=None, me elif product in ['01s', '03s']: resolution = product else: - print (f'ERROR: unknown product {product}. Available only SRTM1 ("01s") and SRTM3 ("03s") DEM using GMT servers') + # open from user-specified NetCDF or GeoTIF file + resolution = None + #print (f'ERROR: unknown product {product}. Available only SRTM1 ("01s") and SRTM3 ("03s") DEM using GMT servers') err, warn = self.validate() #print ('err, warn', err, warn) @@ -207,14 +217,29 @@ def download_dem(self, backend=None, product='SRTM1', resolution_meters=None, me dem_filename = os.path.join(self.basedir, 'DEM_WGS84.nc') with tqdm(desc='DEM Downloading', total=1) as pbar: - ortho = pygmt.datasets.load_earth_relief(resolution=resolution, region=region) + if resolution is not None: + # download DEM + ortho = pygmt.datasets.load_earth_relief(resolution=resolution, region=region) + elif os.path.splitext(product)[-1] in ['.tiff', '.tif', '.TIF']: + ortho = rio.open_rasterio(product, chunks=self.chunksize).squeeze(drop=True)\ + .rename({'y': 'lat', 'x': 'lon'}).sel(lon=slice(minx, maxx))\ + .drop('spatial_ref') + if ortho.lat.diff('lat')[0].item() > 0: + ortho = ortho.sel(lat=slice(miny, maxy)) + else: + ortho = ortho.sel(lat=slice(maxy, miny)) + ortho = ortho.reindex(lat=ortho.lat[::-1]) + elif os.path.splitext(product)[-1] in ['.nc', '.netcdf', '.grd']: + ortho = xr.open_dataarray(product, engine=self.engine, chunks=self.chunksize)\ + .sel(lat=slice(miny, maxy), lon=slice(minx, maxx)) + #print ('ortho', ortho) geoid = xr.open_dataarray(geoid_filename).rename({'y': 'lat', 'x': 'lon'}).interp_like(ortho, method='cubic') #print ('geoid', geoid) if os.path.exists(dem_filename): os.remove(dem_filename) #print ('(ortho + geoid_resamp)', (ortho + geoid)) - (ortho + geoid).rename('dem').to_netcdf(dem_filename) + (ortho + geoid).rename('dem').to_netcdf(dem_filename, engine=self.engine) pbar.update(1) self.dem_filename = dem_filename