#!/usr/bin/env python3
# these are functions for decomposition into E,U vectors from 2 or more tracks
import subprocess as subp
import numpy as np
import xarray as xr
from pathlib import Path
import glob
from scipy import interpolate
from lics_unwrap import *
try:
import dask.array as da
except:
print('warning, no dask installed - the optional function will not work (no worries though)')
import LiCSBAS_inv_lib as inv_lib
import pandas as pd
#from python.get_dates_scihub import startdate
'''
# 2022-10-18 starts here
frame_desc = '082D_05128_030500'
frame_asc = '002A_05136_020502'
asctif='saojorge_offset_A.tif'
desctif='saojorge_offset_D.tif'
beta=25.83
inc_asc, heading_asc = get_frame_inc_heading(frame_asc)
inc_desc, heading_desc = get_frame_inc_heading(frame_desc)
asc = load_tif2xr(asctif)
desc = load_tif2xr(desctif)
# now transform it towards asc:
cube=xr.Dataset()
cube['asc'] = asc
cube['desc'] = desc.interp_like(asc, method='linear'); desc=None
cube['asc_inc'] = inc_asc.interp_like(asc, method='linear'); inc_asc=None
cube['desc_inc'] = inc_desc.interp_like(asc, method='linear'); inc_desc=None
cube['asc_heading'] = heading_asc.interp_like(asc, method='linear'); heading_asc=None
cube['desc_heading'] = heading_desc.interp_like(asc, method='linear'); heading_desc=None
# and decompose:
cube['U']=cube.asc.copy()
cube['E']=cube.asc.copy()
cube['U'].values, cube['E'].values = decompose_np(cube.asc, cube.desc, cube.asc_heading, cube.desc_heading, cube.asc_inc, cube.desc_inc)
'''
[docs]
def calculate_dops_frames(framelist, lon, lat):
''' will calculate dops from a set of frames - now only for given coordinate, can improve if needed
returns: PDOP, HDOP_E, HDOP_N, VDOP '''
from daz_lib import calculate_dops
''' note, for nisar, we can do:
e='unitx.wgs84.tif'
u='unitz.wgs84.tif'
inc, head = extract_inc_heading(e, u, left_looking=True)
but note for azis you need +90!!
'''
# first, get inc and heading:
incs = []
heads = []
for fr in framelist:
inc, head = get_frame_inc_heading(fr)
incs.append(float(inc.sel(lon=lon, lat=lat, method='nearest')))
heads.append(float(head.sel(lon=lon, lat=lat, method='nearest')))
heads = np.array(heads)
incs = np.array(incs)
elevs = 90-incs
azis = heads-90
PDOP, HDOP_E, HDOP_N, VDOP = calculate_dops(elevs, azis)
print("PDOP, HDOP_E, HDOP_N, VDOP = ")
print(PDOP, HDOP_E, HDOP_N, VDOP)
return PDOP, HDOP_E, HDOP_N, VDOP
[docs]
def decompose_framencs(framencs, extract_cum = False, medianfix = False, annual = False,
annual_buffer_months = 0, selperiods = None, do_velUN=False,
do_ENU = False, velname='vel', stdname = None):
""" will decompose frame licsbas results
the basenames in framencs should contain frame id, followed by '.', e.g.:
framencs = ['062D_07629_131313.nc', '172A_07686_131012.nc']
Args:
framencs (list): licsbas nc result files, named by their frame id
extract_cum (bool): if True, will use the first frame and convert to pseudo vertical
annual (bool): if True, will estimate and decompose annual velocities
annual_buffer_months (int): adds extra months for annual velocities
selperiods (list or None): used only if annual=True; override the selection by providing list as [[np.datetime64('2014-01-01'), np.datetime64('2024-01-01')]]
do_velUN ... see decompose_np
do_ENU (bool): if True, it will try invert to ENU rather than fiddling with N (it will also set do_velUN to False)
velname (str): name of the layer to decompose (usually 'vel')
stdname (str): name of the 1-sigma layer to weight velname during decomposition (None means not use)
Returns:
xr.Dataset with U, E, [cum_vert] arrays
"""
interpmethod = 'nearest'
framesetvel = []
frameset = []
firstrun = True
# getting years in all ncs:
yearsall = None
if selperiods:
annual = True # just a workaround..
if do_ENU:
print('using do_ENU - experimental - annuals etc. not implemented yet')
annual=False
selperiods = None
if do_velUN:
print('WARNING - setting do_velUN to False because do_ENU is set to True')
do_velUN = False
for nc in framencs:
framenc = xr.open_dataset(nc)
if 'time' in framenc:
years = framenc.time.dt.year.values
years = list(set(years))
# print(years)
if not yearsall:
yearsall = years
else:
for y in yearsall.copy():
if y not in years:
yearsall.remove(y)
else:
yearsall = [0,1,2]
if len(yearsall)==0:
print('no overlapping year, cancelling annual decomposition')
return False
for nc in framencs:
frame = os.path.basename(nc).split('.')[0]
print('extracting frame '+frame)
inc, heading = get_frame_inc_heading(frame)
framenc = xr.open_dataset(nc)
framevel = framenc[velname]
if medianfix:
framevel = framevel - framevel.median()
if firstrun:
template = framevel.copy()
firstrun = False
if extract_cum:
cum_vert = framenc['cum']
inc = inc.interp_like(framevel, method=interpmethod)
cum_vert = cum_vert/np.cos(np.radians(inc))
else:
framevel = framevel.interp_like(template, method=interpmethod)
if annual:
framenc = framenc.interp_like(template, method=interpmethod)
inc = inc.interp_like(framevel, method=interpmethod)
heading = heading.interp_like(framevel, method=interpmethod)
input_data_set = [framevel.values, heading.values, inc.values]
if stdname:
input_data_set.append(framenc[stdname].interp_like(template, method=interpmethod).values)
framesetvel.append(input_data_set)
if annual:
# doing the annuals!
nc1 = calculate_annual_vels(framenc, yearsall, annual_buffer_months, selperiods)
frameset.append((nc1['vel_annual'], heading.values, inc.values))
dec = xr.Dataset()
U = template.copy()
E = template.copy()
Ustd = template.copy()
Estd = template.copy()
if do_ENU:
N = template.copy()
Nstd = template.copy()
U.values, E.values, N.values, Ustd.values, Estd.values, Nstd.values = decompose_np_multi(framesetvel, do_ENU=do_ENU)
else:
U.values, E.values, Ustd.values, Estd.values = decompose_np_multi(framesetvel, do_velUN=do_velUN)
dec['U'] = U
dec['E'] = E
dec['Ustd'] = Ustd
dec['Estd'] = Estd
if do_ENU:
dec['N'] = N
dec['Nstd'] = Nstd
if annual:
# if annual, then frameset is from nc.vel_annual, heading.values, inc.values
years = None
for framedata in frameset:
if not years:
years = list(framedata[0].year.values)
else:
yearst = list(framedata[0].year.values)
for year in years:
if year not in yearst:
years.remove(year)
# ok, now then decompose the annuals per year
vUxr = frameset[0][0].sel(year = years).copy().rename('vU')*0
vExr = frameset[0][0].sel(year = years).copy().rename('vE')*0
vUstdxr = frameset[0][0].sel(year=years).copy().rename('vUstd') * 0
vEstdxr = frameset[0][0].sel(year=years).copy().rename('vEstd') * 0
for year in years:
annualset = []
for framedata in frameset:
vel = framedata[0].sel(year = year).values
annualset.append((vel, framedata[1], framedata[2]))
print('decomposing year '+str(year))
try:
vU, vE, vUstd, vEstd = decompose_np_multi(annualset, do_velUN=do_velUN) #, beta = 0)
vUxr.loc[year,:,:] = vU
vExr.loc[year,:,:] = vE
vUstdxr.loc[year, :, :] = vUstd
vEstdxr.loc[year, :, :] = vEstd
except:
print('error decomposing, setting nans')
dec['vU'] = vUxr
dec['vE'] = vExr
dec['vUstd'] = vUstdxr
dec['vEstd'] = vEstdxr
if extract_cum:
dec['cum'] = cum_vert
return dec
# Copilot-generated function for custom resample:
# cube is either xr.Dataset or xr.Dataarray
[docs]
def custom_annual_resample(cube, buffermonths=6, buffer_from_midyear = False):
''' Gets annual data resampled per year +- buffermonths
if buffer_from_midyear, it will set buffermonths from the mid-year (June) rather than from Jan and Dec
'''
start_date = pd.to_datetime(cube.time.values[0])
end_date = pd.to_datetime(cube.time.values[-1])
#
# Generate a new time range that extends 6 months before and after
#new_time_range = pd.date_range(
# start=start_date - pd.DateOffset(months=buffermonths),
# end=end_date + pd.DateOffset(months=buffermonths),
# freq='M'
#)
#
# Create an empty list to hold the new resampled cubes
# resampled_cubes = []
# Create an empty dictionary to hold time labels and corresponding values
resampled_data = {'yeardt': [], 'yearvalues': []}
#
for date in pd.date_range(start=start_date, end=end_date, freq='AS'):
# Define the range for the current resample
start_period = date - pd.DateOffset(months=buffermonths)
end_period = date + pd.DateOffset(months=12+buffermonths) - pd.DateOffset(days=1)
if buffer_from_midyear:
start_period = start_period + pd.DateOffset(months=6)
end_period = end_period - pd.DateOffset(months=6)
#
# Select the data within this range
selected_data = cube.sel(time=slice(start_period, end_period))
#
# Create a new time coordinate for the selected data period
new_time = pd.date_range(start=start_period, end=end_period, freq='M')
selected_data = selected_data.reindex({'time': new_time}, method='nearest')
#
# Append the selected data to the list
#resampled_cubes.append(selected_data)
resampled_data['yeardt'].append(date)
resampled_data['yearvalues'].append(selected_data)
#
# Combine all resampled data into one DataArray or Dataset
#resampled_cube = xr.concat(resampled_cubes, dim='time')
#return resampled_cube
#
# Convert lists to pandas DataFrame or xarray Dataset for easier use
yeardt = pd.to_datetime(resampled_data['yeardt'])
yearvalues = xr.concat(resampled_data['yearvalues'], dim=pd.Index(yeardt, name='time'))
#
return yeardt, yearvalues
[docs]
def calculate_annual_vels(cube, commonyears = None, buffermonths = 0, selperiods = None):
"""Will calculate annual velocities from LiCSBAS results
Args:
cube (xr.Dataset): loaded netcdf file, extracted using e.g. LiCSBAS_out2nc.py
commonyears (list): list of years to decompose
buffermonths (int): extend selection of annual data by a number of +-buffermonths months (experimental)
selperiods (list or None): override the selection by providing list as [[np.datetime64('2014-01-01'), np.datetime64('2024-01-01')]]
Returns:
xr.Dataset with new dataarray: vel_annual
"""
if commonyears:
cube = cube.sel(time=np.isin(cube.time.dt.year.values, commonyears))
if type(selperiods) != type(None):
annualset = []
for sp in selperiods:
startdate = sp[0]
enddate = sp[1]
yeardt = startdate+(enddate-startdate)/2+5
yearcum = cube['cum'].sel(time=slice(startdate, enddate))
annualset.append([yeardt, yearcum])
elif buffermonths > 0:
print('Warning, using Copilot-generated trick to add more months around year of interest...')
annualset = custom_annual_resample(cube['cum'], buffermonths)
else:
annualset = cube.cum.resample(time='AS')
firstrun = True
for yeardt, yearcum in annualset:
year = str(yeardt).split('-')[0]
print('processing year '+year)
dt_cum = (np.array([tdate.toordinal() for tdate in yearcum.time.dt.date.values]) - pd.Timestamp(yeardt).date().toordinal())/365.25 # in fraction of a year
# see LiCSBAS_cum2vel.py:
cum_tmp = yearcum.values
n_im, length, width = cum_tmp.shape
bool_allnan = np.all(np.isnan(cum_tmp), axis=0)
vconst = np.zeros((length, width), dtype=np.float32)*np.nan
vel = np.zeros((length, width), dtype=np.float32)*np.nan
#
cum_tmp = cum_tmp.reshape(n_im, length*width)[:, ~bool_allnan.ravel()].transpose()
vel[~bool_allnan], vconst[~bool_allnan] = inv_lib.calc_vel(cum_tmp, dt_cum)
#vel[~bool_allnan], vconst[~bool_allnan] = calc_vel(cum_tmp, dt_cum)
vel_annual = cube.vel.copy()
vel_annual.values = vel
vel_annual = vel_annual.assign_coords({'year':int(year)}).expand_dims('year').rename('vel_annual')
if firstrun:
vel_annual_cube = vel_annual.copy()
firstrun = False
#dv = vel
else:
vel_annual_cube = xr.concat([vel_annual_cube,vel_annual], dim='year')
#dv = np.vstack(dv,vel)
cube['vel_annual'] = vel_annual_cube.copy()
return cube
[docs]
def get_frame_inc_heading(frame):
"""Extracts inc and heading from E, U for given frame
"""
geoframedir = os.path.join(os.environ['LiCSAR_public'], str(int(frame[:3])), frame)
# look angle (inc) / heading - probably ok, but needs check:
e=os.path.join(geoframedir,'metadata',frame+'.geo.E.tif')
#n=os.path.join(geoframedir,'metadata',frame+'.geo.N.tif') #no need for N
u=os.path.join(geoframedir,'metadata',frame+'.geo.U.tif')
return extract_inc_heading(e, u)
def get_frame_enu(frame, template = None):
geoframedir = os.path.join(os.environ['LiCSAR_public'], str(int(frame[:3])), frame)
# look angle (inc) / heading - probably ok, but needs check:
e = os.path.join(geoframedir, 'metadata', frame + '.geo.E.tif')
n = os.path.join(geoframedir,'metadata',frame+'.geo.N.tif') #no need for N ?
u = os.path.join(geoframedir, 'metadata', frame + '.geo.U.tif')
u = load_tif2xr(u)
e = load_tif2xr(e)
n = load_tif2xr(n)
if type(template) != type(None):
u = u.interp_like(template, method='nearest')
e = e.interp_like(template, method='nearest')
n = n.interp_like(template, method='nearest')
if np.isnan(u.mean()):
u = e * 0 # just in case... for bovls.. although probably not needed
return e, n, u
[docs]
def decompose_dask(cube, blocklen=5, num_workers=5):
"""Simple parallel decomposition of dec. datacube (must have asc,desc,asc_inc,desc_inc, asc_heading, desc_heading arrays)
"""
winsize = (blocklen, blocklen)
asc = da.from_array(cube['asc'].astype(np.float32), chunks=winsize)
desc = da.from_array(cube['desc'].astype(np.float32), chunks=winsize)
ascinc = da.from_array(cube['asc_inc'].astype(np.float32), chunks=winsize)
descinc = da.from_array(cube['desc_inc'].astype(np.float32), chunks=winsize)
aschead = da.from_array(cube['asc_heading'].astype(np.float32), chunks=winsize)
deschead = da.from_array(cube['desc_heading'].astype(np.float32), chunks=winsize)
#f = da.map_blocks(decompose_np, asc, desc, aschead, deschead, ascinc, descinc, beta=0, meta=np.array((),())) #, chunks = (1,1))
f = da.map_blocks(decompose_np, asc, desc, aschead, deschead, ascinc, descinc, beta=0, meta=(np.array((), dtype=np.float32), np.array((), dtype=np.float32)))
return f.compute(num_workers=num_workers)
[docs]
def decompose_xr(asc, desc, heading_asc, heading_desc, inc_asc, inc_desc, do_velUN = True):
"""Perform simple decomposition for two frames in asc and desc.
inputs are xr.dataarrays - this will also check/interpolate them to fit
Note - better use decompose_np_multi to get also sigmas etc.
"""
cube=xr.Dataset()
cube['asc'] = asc
cube['desc'] = desc.interp_like(asc, method='nearest'); desc=None
cube['U']=cube.asc.copy()
cube['E']=cube.asc.copy()
if not np.isscalar(heading_asc):
cube['asc_heading'] = heading_asc.interp_like(asc, method='linear'); heading_asc=cube.asc_heading.values
cube['desc_heading'] = heading_desc.interp_like(asc, method='linear'); heading_desc=cube.desc_heading.values
if not np.isscalar(inc_asc):
cube['asc_inc'] = inc_asc.interp_like(asc, method='linear'); inc_asc=cube.asc_inc.values
cube['desc_inc'] = inc_desc.interp_like(asc, method='linear'); inc_desc=cube.desc_inc.values
cube['U'].values, cube['E'].values = decompose_np(cube.asc.values, cube.desc.values, heading_asc, heading_desc, inc_asc , inc_desc, do_velUN = do_velUN)
return cube[['U', 'E']]
# 2022-10-18 - this should be pretty good one (next only use weights or something)
[docs]
def decompose_np(vel_asc, vel_desc, aschead, deschead, ascinc, descinc, beta=0, do_velUN = True):
"""Decomposes values from ascending and descending np (or xr) arrays, using heading and inc. angle
(these might be arrays of same size of just float values)
Args:
beta (float): angle of expected horizontal motion direction, clockwise from the E, in degrees
do_velUN (boolean): if yes, extract v_{UN} instead of v_U, following Qi et al, 2022: doi=10.1029/2F2022JB024176
"""
vel_E = np.zeros(vel_desc.shape)
vel_U = np.zeros(vel_desc.shape)
#
if do_velUN:
U_asc = np.sqrt(1 - (np.sin(np.radians(ascinc))**2) * (np.cos(np.radians(aschead))**2))
U_desc = np.sqrt(1 - (np.sin(np.radians(descinc)) ** 2) * (np.cos(np.radians(deschead)) ** 2))
else:
U_asc = np.cos(np.radians(ascinc))
U_desc = np.cos(np.radians(descinc))
#
E_asc = -np.sin(np.radians(ascinc))*np.cos(np.radians(aschead+beta))
E_desc = -np.sin(np.radians(descinc))*np.cos(np.radians(deschead+beta))
#
for ii in np.arange(0,vel_E.shape[0]):
for jj in np.arange(0,vel_E.shape[1]):
# velocities
d = np.array([[vel_asc[ii,jj], vel_desc[ii,jj]]]).T
# if the velocities contain nan, will return nan:
if np.isnan(np.max(d)):
vel_U[ii,jj] = np.nan
vel_E[ii,jj] = np.nan
else:
# create the design matrix
if np.isscalar(U_asc): # in case of only values (i.e. one inc and heading per each frame)
G = np.array([[U_asc, E_asc], [U_desc, E_desc]])
else: # in case this is array
G = np.array([[U_asc[ii,jj], E_asc[ii,jj]], [U_desc[ii,jj], E_desc[ii,jj]]])
# solve the linear system for the Up and East velocities
m = np.linalg.solve(G, d)
# save to arrays
vel_U[ii,jj] = m[0]
vel_E[ii,jj] = m[1]
return vel_U, vel_E
'''
this is to load 3 datasets and decompose them:
dirpath='/gws/ssde/j25a/nceo_geohazards/vol1/public/shared/temp/earmla'
#for frame in []
nc1 = os.path.join(dirpath, '051D_03973_131313.nc')
nc1=xr.open_dataset(nc1)
vel1 = nc1.vel.values
heading1 = -169.87
inc1 = 43.64
nc2 = os.path.join(dirpath, '124D_04017_131313.nc')
nc2=xr.open_dataset(nc2)
vel2 = nc2.vel.interp_like(nc1.vel).values
heading2 = -169.88
inc2 = 34.98
nc3 = os.path.join(dirpath, '175A_03997_131313.nc')
nc3=xr.open_dataset(nc3)
vel3 = nc3.vel.interp_like(nc1.vel).values
heading3 = -10.16
inc3 = 38.42
years = np.array([2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022])
vUxr = nc1.vel_annual.sel(year = years).copy().rename('vU')
vExr = nc1.vel_annual.sel(year = years).copy().rename('vE')
for year in years:
vel1 = nc1.vel_annual.sel(year = year).values
vel2 = nc2.vel_annual.sel(year = year).values
vel3 = nc3.vel_annual.sel(year = year).values
input_data = [(vel1, heading1, inc1), (vel2, heading2, inc2), (vel3, heading3, inc3)]
print('decomposing year '+str(year))
vU, vE = decompose_np_multi(input_data, beta = 0)
vUxr.loc[year,:,:] = vU
vExr.loc[year,:,:] = vE
decomposedxr = xr.Dataset()
decomposedxr['vU'] = vUxr
decomposedxr['vE'] = vExr
decomposedxr.to_netcdf('decomposed_s1.nc')
'''
[docs]
def decompose_geotiffs(veltifs, Etifs, Utifs, Ntifs = None, vstdtifs = None, leftlooking = None, # aziflags = None,
do_ENU = False, do_velUN = False):
""" Decompose set of geotiffs
inputs are list of respective geotiffs.
vstdtifs and leftlooking can be None to skip.
Similarly, N tifs are not required - but they are necessary if bovls are used! In case of not providing Ns, it will use old approach to estimate from heading/inc. angle
if used, leftlooking must be list, e.g. [False, False, True] meaning the third set is NISAR. Note, leftlooking is used only for method with only E, U to decompose..
Note: aziflags do not work in this scope, so they are skipped:
if bovls are included, you need to provide aziflags such as [None, 'D', 'A', None] meaning second tifs are azimuth (bovls) of descending track
(full flags are also accepted, i.e. ['D','D','A','A'] would work ok even if first and last are not azimuth)
"""
input_data = [] #(vel1, heading1, inc1), (vel2, heading2, inc2), (vel3, heading3, inc3)]
firstrun = True
use_enu_vectors = False
if Ntifs:
use_enu_vectors = True
for i in range(len(veltifs)):
vel = load_tif2xr(veltifs[i])
if firstrun:
template = vel.copy()
firstrun = False
else:
vel = vel.interp_like(template, method='nearest')
if leftlooking:
left = leftlooking[i]
else:
left = False
#if aziflags:
# azi = aziflags[i]
#else:
# azi = None
if not Ntifs:
inc, head = extract_inc_heading(Etifs[i], Utifs[i], left_looking=left) #, aziflag=azi)
inc = inc.interp_like(template, method='nearest')
head = head.interp_like(template, method='nearest')
if vstdtifs:
vstd = load_tif2xr(vstdtifs[i])
vstd = vstd.interp_like(template, method='nearest')
input_data.append((vel.values, head.values, inc.values, vstd.values))
else:
input_data.append((vel.values, head.values, inc.values))
else:
u = load_tif2xr(Utifs[i])
e = load_tif2xr(Etifs[i])
n = load_tif2xr(Ntifs[i])
u = u.interp_like(template, method='nearest')
e = e.interp_like(template, method='nearest')
n = n.interp_like(template, method='nearest')
if np.isnan(u.mean()):
u = e * 0 # just in case... although probably not needed
if vstdtifs:
vstd = load_tif2xr(vstdtifs[i])
vstd = vstd.interp_like(template, method='nearest')
input_data.append((vel.values, e.values, n.values, u.values, vstd.values))
else:
input_data.append((vel.values, e.values, n.values, u.values))
#
velouts = decompose_np_multi(input_data, do_velUN=do_velUN, do_ENU = do_ENU, input_is_enu_vectors = use_enu_vectors)
if do_ENU:
vel_U, vel_E, vel_N, vel_Ustd, vel_Estd, vel_Nstd = velouts
else:
vel_U, vel_E, vel_Ustd, vel_Estd = velouts
# get it back as netcdf:
dec = xr.Dataset()
dec['U'] = template.copy()
dec['E'] = template.copy()
if vstdtifs:
dec['Ustd'] = template.copy()
dec['Estd'] = template.copy()
if do_ENU:
dec['N'] = template.copy()
if vstdtifs:
dec['Nstd'] = template.copy()
#
dec['U'].values = vel_U
dec['E'].values = vel_E
if vstdtifs:
dec['Ustd'].values = vel_Ustd
dec['Estd'].values = vel_Estd
if do_ENU:
dec['N'].values = vel_N
if vstdtifs:
dec['Nstd'].values = vel_Nstd
return dec
[docs]
def decompose_np_multi(input_data, beta = 0, do_velUN=False, do_ENU = False, input_is_enu_vectors = False):
"""Decompose 2 or more frames
Args:
input data (list of tuples) e.g. input_data = [(vel1, heading1, inc1), (vel2, heading2, inc2), (vel3, heading3, inc3)]
... in case there are 4, we will assume vstd (1-sigma) as:
[(vel1, heading1, inc1, vstd1), (vel2, heading2, inc2, vstd2), ...]
do_velUN and do_ENU: see decompose_framencs
do_ENU: returns U, E, N (and Ustd, Estd, Nstd)
input_is_enu_vectors (bool): if True, the input_data is actually expected to be in the form of [(vel1, E1, N1, U1), ...] .. or with vstd, sure
Returns:
4x np.ndarray of decomposed outputs U, E, and their 1-sigma Ustd, Estd
Note: velX is np.array and headingX/incX is in degrees, either a number or np.array
"""
#
template = input_data[0][0]
vel_E = np.zeros(template.shape)
vel_U = np.zeros(template.shape)
vel_Estd = np.zeros(template.shape)
vel_Ustd = np.zeros(template.shape)
#
Us=list()
Es=list()
if do_ENU:
vel_N = np.zeros(template.shape)
vel_Nstd = np.zeros(template.shape)
Ns=list()
do_velUN=False
if input_is_enu_vectors and do_velUN:
print('ERROR: velUN approach starts from heading and inc data to recalc UN. While simple to rearrange, simpler for you to just... use extract_inc_heading or get_frame_inc_heading')
return False
vels = []
vstds = []
if do_velUN and input_is_enu_vectors:
print('WARNING - this setting works ok only if all your inputs are from right-looking satellite (e.g. S1, but not NISAR)')
for frame in input_data:
vel = frame[0]
if input_is_enu_vectors:
E = frame[1]
N = frame[2]
U = frame[3]
if len(frame) > 4:
vstd = frame[4]
else:
vstd = vel * 0 + 1
if do_velUN:
incangle, heading = extract_inc_heading(E, U, eu_are_files=False)
U = np.sqrt(1 - (np.sin(np.radians(incangle)) ** 2) * (np.cos(np.radians(incangle)) ** 2))
Us.append(U)
Es.append(E)
if do_ENU:
Ns.append(N)
vels.append(vel)
vstds.append(vstd)
else:
heading = frame[1]
incangle = frame[2]
if len(frame)>3:
vstd = frame[3]
else:
vstd = vel*0+1
if do_velUN:
U = np.sqrt(1 - (np.sin(np.radians(incangle)) ** 2) * (np.cos(np.radians(incangle)) ** 2))
else:
U = np.cos(np.radians(incangle))
#Us = np.append(Us, np.cos(np.radians(incangle)))
Us.append(U)
#Es = np.append(Es, -np.sin(np.radians(incangle))*np.cos(np.radians(heading+beta)))
Es.append(-np.sin(np.radians(incangle))*np.cos(np.radians(heading+beta)))
vels.append(vel)
vstds.append(vstd)
if do_ENU:
#
Ns.append(np.sin(np.radians(incangle))*np.sin(np.radians(heading+beta)))
# run for each pixel
numframes = len(vels)
Us = np.array(Us)
Es = np.array(Es)
if do_ENU:
Ns = np.array(Ns)
for ii in np.arange(0,vel_E.shape[0]):
for jj in np.arange(0,vel_E.shape[1]):
# prepare template for d = G m
d = np.array(())
Qd = np.array(())
for i in range(numframes):
d = np.append(d, np.array([vels[i][ii,jj]]))
Qd = np.append(Qd, np.array([vstds[i][ii, jj]**2])) # should add variances to Qd
d = np.array([d]).T
Qd=np.diag(Qd)
if np.isnan(d).any():
# if at least one is nan, skip it: # can improve it but 'all' is not an option
#if np.isnan(np.max(d)):
vel_U[ii,jj] = np.nan
vel_E[ii,jj] = np.nan
if do_ENU:
vel_N[ii, jj] = np.nan
else:
# create the design matrix
if np.isscalar(Us[0]): # in case of only values (i.e. one inc and heading per each frame)
if do_ENU:
G = np.vstack([Us, Es, Ns]).T
else:
G = np.vstack([Us, Es]).T
else: # in case this is array # not tested!
if do_ENU:
G = np.vstack([Us[:,ii,jj], Es[:,ii,jj], Ns[:,ii,jj]]).T
else:
G = np.vstack([Us[:,ii,jj], Es[:,ii,jj]]).T
# solve the linear system for the Up and East velocities
#m = np.linalg.solve(G, d)
try:
# 2025/09: thanks A. Watson on his https://github.com/andwatson/decompose_insar_velocities
Qd_inv = np.linalg.inv(Qd) # this means weights..
# m
m = np.linalg.inv(np.dot(G.T, np.dot(Qd_inv, G))) @ np.dot(G.T, np.dot(Qd_inv, d))
# Qm
Qm = np.linalg.inv(np.dot(G.T, np.dot(Qd_inv, G)))
# The lstsq solution does not give Qm ...
#m = np.linalg.lstsq(G, d / np.sqrt(Qd))[0]
#m = m*np.sqrt(Qd)
#m = m[:,0]
# save to arrays
vel_U[ii,jj] = m[0].item()
vel_E[ii,jj] = m[1].item()
vel_Ustd[ii,jj] = np.sqrt(Qm[0,0])
vel_Estd[ii, jj] = np.sqrt(Qm[1,1])
if do_ENU:
vel_N[ii,jj] = m[2].item()
vel_Nstd[ii,jj] = np.sqrt(Qm[2,2])
except ValueError as e:
print(f"Error: {e}. Setting nan")
vel_U[ii,jj] = np.nan
vel_E[ii,jj] = np.nan
vel_Ustd[ii, jj] = np.nan
vel_Estd[ii, jj] = np.nan
if do_ENU:
vel_N[ii,jj] = np.nan
vel_Nstd[ii, jj] = np.nan
vel_Ustd[vel_Ustd == 0] = np.nan
vel_Estd[vel_Estd == 0] = np.nan
if do_ENU:
vel_Nstd[vel_Nstd == 0] = np.nan
return vel_U, vel_E, vel_N, vel_Ustd, vel_Estd, vel_Nstd
else:
return vel_U, vel_E, vel_Ustd, vel_Estd
'''
files = [
"022D.nc",
"146A.nc",
"095D.nc",
]
# or
import glob
# files = glob.glob('*.nc')
'''
[docs]
def decompose_cum_framencs(framencs, time_step_days=14, do_velUN=False, do_ENU=False):
''' This will decompose cum values into given time steps.
Args:
framencs (list): list of input nc files - note, naming must be the frame IDs, e.g. ['022D_03989_131313.nc', '044A_03932_131313.nc', ...]
time_step_days (int): output time step of decomposed components
do_velUN (bool): if True, it will merge U+N into returned 'pseudovertical' - use only with standard LOS data from right-looking sats (i.e. not NISAR-ready)
do_ENU (bool): if True, will perform full decomposition
Returns:
xr.Dataset
'''
if do_velUN and do_ENU:
print('ERROR - please use only one of the options')
return False
ds_out = regrid_netcdf_collection(
framencs,
var="cum",
time_step_days=time_step_days,
time_method="nearest",
time_tolerance_days=21,
)
print(ds_out)
# ds_out.to_netcdf('regridded_cum.nc')
# now, load ENU files:
template = ds_out['cum'][0][0]
enu = {
frame: get_frame_enu(frame, template)
for frame in ds_out.frame.values
}
ds_out["E"] = xr.concat(
[enu[f][0].expand_dims(frame=[f]) for f in ds_out.frame.values],
dim="frame",
)
ds_out["N"] = xr.concat(
[enu[f][1].expand_dims(frame=[f]) for f in ds_out.frame.values],
dim="frame",
)
ds_out["U"] = xr.concat(
[enu[f][2].expand_dims(frame=[f]) for f in ds_out.frame.values],
dim="frame",
)
nt = ds_out.sizes["time"]
ny = ds_out.sizes["lat"]
nx = ds_out.sizes["lon"]
cum_U = np.full((nt, ny, nx), np.nan, dtype=np.float32)
cum_E = np.full((nt, ny, nx), np.nan, dtype=np.float32)
cum_U_std = np.full((nt, ny, nx), np.nan, dtype=np.float32)
cum_E_std = np.full((nt, ny, nx), np.nan, dtype=np.float32)
if do_ENU:
cum_N = np.full((nt, ny, nx), np.nan, dtype=np.float32)
cum_N_std = np.full((nt, ny, nx), np.nan, dtype=np.float32)
E = ds_out.E.values
N = ds_out.N.values
U = ds_out.U.values
if do_velUN:
inc, head = zip(*[
extract_inc_heading(
E[i], U[i],
left_looking=False,
eu_are_files=False
)
for i in range(E.shape[0])
])
inc = np.stack(inc)
head = np.stack(head)
for t in range(nt):
print('Decomposing '+str(t+1)+'/'+str(nt))
cum_t = ds_out.cum.isel(time=t).values # (frame, lat, lon)
if do_velUN:
input_data = [
(cum_t[f], head[f], inc[f])
for f in range(cum_t.shape[0])
]
else:
input_data = [
(cum_t[f], E[f], N[f], U[f])
for f in range(cum_t.shape[0])
]
outdec = decompose_np_multi(input_data,
do_velUN=do_velUN,
do_ENU=do_ENU, # ... this will be great...
input_is_enu_vectors = not do_velUN)
if do_ENU:
cum_U[t], cum_E[t], cum_N[t], cum_U_std[t], cum_E_std[t], cum_N_std[t] = outdec
else:
cum_U[t], cum_E[t], cum_U_std[t], cum_E_std[t] = outdec
coords = {
"time": ds_out.time,
"lat": ds_out.lat,
"lon": ds_out.lon,
}
if do_velUN:
outcumu = 'cum_UN'
else:
outcumu = 'cum_U'
ds_out[outcumu] = xr.DataArray(
cum_U,
coords=coords,
dims=("time", "lat", "lon"),
)
ds_out["cum_E"] = xr.DataArray(
cum_E,
coords=coords,
dims=("time", "lat", "lon"),
)
ds_out[outcumu+"_std"] = xr.DataArray(
cum_U_std,
coords=coords,
dims=("time", "lat", "lon"),
)
ds_out["cum_E_std"] = xr.DataArray(
cum_E_std,
coords=coords,
dims=("time", "lat", "lon"),
)
if do_ENU:
ds_out["cum_N"] = xr.DataArray(
cum_E,
coords=coords,
dims=("time", "lat", "lon"),
)
ds_out["cum_N_std"] = xr.DataArray(
cum_U_std,
coords=coords,
dims=("time", "lat", "lon"),
)
return ds_out
[docs]
def regrid_netcdf_collection(
files,
var="cum",
time_step_days=14,
spatial_ref=None,
time_method="nearest",
time_tolerance_days=None,
):
"""
Parameters
----------
files : list[str]
List of NetCDF files.
var : str
Data variable name.
time_step_days : int
Output timestep.
spatial_ref : str | None
File defining output lat/lon grid.
If None, first file is used.
time_method : str
'nearest' or 'linear'.
time_tolerance_days : float | None
Maximum temporal distance allowed for nearest interpolation.
Returns
-------
xr.Dataset
Dataset containing:
cum(frame,time,lat,lon)
and coordinate:
frame = ['022D','146A',...]
"""
# ----------------------------
# open datasets
# ----------------------------
dsets = [xr.open_dataset(f) for f in files]
frame_names = [
Path(f).stem
for f in files
]
# ----------------------------
# common temporal overlap
# ----------------------------
overlap_start = max(
pd.Timestamp(ds.time.min().values)
for ds in dsets
)
overlap_end = min(
pd.Timestamp(ds.time.max().values)
for ds in dsets
)
if overlap_end <= overlap_start:
raise ValueError(
"No common temporal overlap between datasets."
)
common_time = pd.date_range(
overlap_start,
overlap_end,
freq=f"{time_step_days}D"
)
print('Output min dt: '+str(overlap_start))
print('Output max dt: ' + str(overlap_end))
# ----------------------------
# reference spatial grid
# ----------------------------
if spatial_ref is None:
ref_ds = dsets[0]
else:
ref_ds = xr.open_dataset(spatial_ref)
target_grid = xr.Dataset(
coords={
"lat": ref_ds.lat,
"lon": ref_ds.lon,
}
)
frames = []
for ds, frame in zip(dsets, frame_names):
print('Processing frame '+frame)
da = ds[var]
# ---------------------------------
# spatial nearest-neighbour regrid
# ---------------------------------
print('...regridding in space')
da = da.interp(
lat=target_grid.lat,
lon=target_grid.lon,
method="nearest",
)
# ---------------------------------
# temporal interpolation
# ---------------------------------
# remove completely-empty timesteps
valid = ~da.isnull().all(dim=("lat", "lon"))
da = da.sel(time=valid)
print('...interpolating in time - using method '+time_method)
if time_method == "nearest":
if time_tolerance_days is None:
da = da.reindex(
time=common_time,
method="nearest",
)
else:
da = da.reindex(
time=common_time,
method="nearest",
tolerance=np.timedelta64(
time_tolerance_days, "D"
),
)
elif time_method == "linear":
da = da.interp(
time=common_time,
method="linear",
)
else:
raise ValueError(
"time_method must be 'nearest' or 'linear'"
)
da = da.expand_dims(frame=[frame])
frames.append(da)
out = xr.concat(frames, dim="frame")
return xr.Dataset(
{
var: out
}
)
''' plotting decomposed cum:
import numpy as np
import matplotlib.pyplot as plt
lat0, lon0 = 50.127958, 14.020140
ts = ds_out["cum_U"].sel(lat=lat0, lon=lon0, method="nearest")
ts_std = ds_out["cum_U_std"].sel(lat=lat0, lon=lon0, method="nearest")
# or
lat_min, lat_max = 50.123377, 50.131593
lon_min, lon_max = 14.010899, 14.026722
roi = ds_out.sel(
lat=slice(lat_min, lat_max),
lon=slice(lon_min, lon_max)
)
ts = roi.cum_U.mean(("lat", "lon"), skipna=True)
ts_std = roi.cum_U.std(
("lat", "lon"),
skipna=True
)
t = (
(ts.time.values - ts.time.values[0])
/ np.timedelta64(1, "D")
) / 365.25
G = np.column_stack([
np.ones_like(t),
t,
np.sin(2*np.pi*t),
np.cos(2*np.pi*t),
])
W = np.diag(1/ts_std.values**2)
mask = np.isfinite(ts.values) & np.isfinite(ts_std.values)
Gm = G[mask]
ym = ts.values[mask]
Wm = np.diag(1/ts_std.values[mask]**2)
m = np.linalg.solve(
Gm.T @ Wm @ Gm,
Gm.T @ Wm @ ym
)
model = G @ m
plt.figure(figsize=(10,5))
plt.plot(ts.time, ts, "ko", ms=4, label="cum_U")
plt.fill_between(
ts.time.values,
ts.values - ts_std.values,
ts.values + ts_std.values,
alpha=0.3,
color="C0",
label="±1σ",
)
plt.plot(ts.time, model, "r-", lw=2, label="model")
plt.ylim(-5,5)
plt.ylabel('mm')
plt.title(str(lon0)+', '+str(lat0))
plt.legend()
plt.grid()
plt.show()
'''
'''
aschead=-9.918319
deschead=-169.61931
ascinc=39.4824
descinc=33.7491
desc=xr.open_dataset('082D_05128_030500ok.nc')
asc=xr.open_dataset('002A_05136_020502.nc')
D=desc.cum[-2]-desc.cum[-3]
A=asc.cum[-1]-asc.cum[-4]
from unwrp_multiscale import *
export_xr2tif(A,'A.tif', dogdal=False)
export_xr2tif(D,'D.tif', dogdal=False)
os.system('gdalwarp2match.py A.tif D.tif Aok.tif')
os.system('gdalwarp2match.py D.tif Aok.tif Dok.tif')
'''
''' (old) usage example:
#def decompose_xr(asc, desc, aschead, deschead, ascinc, descinc):
#
U,E = decompose('Aok.tif', 'Dok.tif', aschead, deschead, ascinc, descinc)
aa = rioxarray.open_rasterio('Aok.tif')
aa.values[0]=U
export_xr2tif(aa,'U.tif', lonlat=False,dogdal=False)
import LiCSBAS_io_lib as io
# to load to np
asctif=...
desctif=...
vel_asc = io.read_geotiff(asctif)
vel_desc = io.read_geotiff(desctif)
aschead=349.613
deschead=190.3898
ascinc=34.403
descinc=34.240
'''
'''
orig AW approach:
import matplotlib.pyplot as plt
# these packages are only needed for the final multivariate plot
import seaborn as sns
import pandas as pd
import interseis_lib as lib
# setup file names
vel_file_asc = 'data/087A_04904_121313_vel'
par_file_asc = 'data/087A_04904_121313.par'
E_file_asc = 'data/087A_04904_121313_E.geo'
N_file_asc = 'data/087A_04904_121313_N.geo'
U_file_asc = 'data/087A_04904_121313_U.geo'
vel_file_desc = 'data/167D_04884_131212_vel'
par_file_desc = 'data/167D_04884_131212.par'
E_file_desc = 'data/167D_04884_131212_E.geo'
N_file_desc = 'data/167D_04884_131212_N.geo'
U_file_desc = 'data/167D_04884_131212_U.geo'
# read array dimensions from par file
width_asc = int(lib.get_par(par_file_asc,'width'))
length_asc = int(lib.get_par(par_file_asc,'nlines'))
width_desc = int(lib.get_par(par_file_desc,'width'))
length_desc = int(lib.get_par(par_file_desc,'nlines'))
# get corner positions
corner_lat_asc = float(lib.get_par(par_file_asc,'corner_lat'))
corner_lon_asc = float(lib.get_par(par_file_asc,'corner_lon'))
corner_lat_desc = float(lib.get_par(par_file_desc,'corner_lat'))
corner_lon_desc = float(lib.get_par(par_file_desc,'corner_lon'))
# get post spacing (distance between velocity measurements)
post_lat_asc = float(lib.get_par(par_file_asc,'post_lat'))
post_lon_asc = float(lib.get_par(par_file_asc,'post_lon'))
post_lat_desc = float(lib.get_par(par_file_desc,'post_lat'))
post_lon_desc = float(lib.get_par(par_file_desc,'post_lon'))
# calculate grid spacings
lat_asc = corner_lat_asc + post_lat_asc*np.arange(1,length_asc+1) - post_lat_asc/2
lon_asc = corner_lon_asc + post_lon_asc*np.arange(1,width_asc+1) - post_lon_asc/2
lat_desc = corner_lat_desc + post_lat_desc*np.arange(1,length_desc+1) - post_lat_desc/2
lon_desc = corner_lon_desc + post_lon_desc*np.arange(1,width_desc+1) - post_lon_desc/2
# load in velocities
vel_asc = np.fromfile(vel_file_asc, dtype='float32').reshape((length_asc, width_asc))
vel_desc = np.fromfile(vel_file_desc, dtype='float32').reshape((length_desc, width_desc))
# load in unit vectors
E_asc = np.fromfile(E_file_asc, dtype='float32').reshape((length_asc, width_asc))
N_asc = np.fromfile(N_file_asc, dtype='float32').reshape((length_asc, width_asc))
U_asc = np.fromfile(U_file_asc, dtype='float32').reshape((length_asc, width_asc))
E_desc = np.fromfile(E_file_desc, dtype='float32').reshape((length_desc, width_desc))
N_desc = np.fromfile(N_file_desc, dtype='float32').reshape((length_desc, width_desc))
U_desc = np.fromfile(U_file_desc, dtype='float32').reshape((length_desc, width_desc))
# load the naf fault trace
fault_trace = np.loadtxt('data/naf_trace.xy')
# pre-allocate
vel_E = np.zeros((len(lat_regrid), len(lon_regrid)))
vel_U = np.zeros((len(lat_regrid), len(lon_regrid)))
# loop through every pixel
for ii in np.arange(0,len(lat_regrid)):
for jj in np.arange(0,len(lon_regrid)):
# create the design matrix
G = np.array([[U_asc_regrid[ii,jj], E_asc_regrid[ii,jj]], [U_desc_regrid[ii,jj], E_desc_regrid[ii,jj]]])
# get the two velocities for this pixel
d = np.array([[vel_asc_regrid[ii,jj], vel_desc_regrid[ii,jj]]]).T
# solve the linear system for the Up and East velocities
m = np.linalg.solve(G, d)
# save to arrays
vel_U[ii,jj] = m[0]
vel_E[ii,jj] = m[1]
'''