Skip to content

Commit

Permalink
Allow apply flat-earth correction without topo phase using topo=0 (to…
Browse files Browse the repository at this point in the history
…po=None means no corrections).
  • Loading branch information
AlexeyPechnikov committed Feb 28, 2024
1 parent e72441f commit 982a770
Showing 1 changed file with 5 additions and 2 deletions.
7 changes: 5 additions & 2 deletions pygmtsar/pygmtsar/Stack_phasediff.py
Original file line number Diff line number Diff line change
Expand Up @@ -496,15 +496,18 @@ def block_phasediff(date1, date2, prm1, prm2, ylim, xlim):

# use outer variables topo, data1, data2, prm1, prm2
# build topo block
dy, dx = topo.y.diff('y').item(0), topo.x.diff('x').item(0)
if dy == 1 and dx == 1:
if not isinstance(topo, xr.DataArray):
# topography is a constant, typically, zero
block_topo = topo * xr.ones_like(block_data1, dtype=np.float32)
elif dy == 1 and dx == 1:
# topography is already in the original resolution
block_topo = topo.isel(y=slice(ylim[0], ylim[1]), x=slice(xlim[0], xlim[1]))\
.compute(n_workers=1)\
.fillna(0)\
.assign_coords(y=ys, x=xs)
else:
# topography resolution is different, interpolation with extrapolation required
dy, dx = topo.y.diff('y').item(0), topo.x.diff('x').item(0)
# convert indices 0.5, 1.5,... to 0,1,... for easy calculations
# fill NaNs by zero because typically DEM is missed outside of land areas
block_topo = topo.sel(y=slice(ys[0]-2*dy, ys[-1]+2*dy), x=slice(xs[0]-2*dx, xs[-1]+2*dx))\
Expand Down

0 comments on commit 982a770

Please sign in to comment.