Yang has provide this code for patching
#!/usr/bin/env python3
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Yang Lei, Jet Propulsion Laboratory
#November 2017
import argparse
import isce
import isceobj
import numpy as np
import shelve
import os
import datetime
from isceobj.Constants import SPEED_OF_LIGHT
import pdb
import logging
from imageMath import IML
from scipy import signal
import scipy.io as sio
import rioxarray
import xarray
from osgeo import osr,gdal
def cmdLineParse():
'''
Command line parser.
'''
parser = argparse.ArgumentParser( description='Converting netcdf variables to geotiff')
parser.add_argument('-i', '--nc', dest='nc', type=str, required=True,
help = 'input netcdf file path')
parser.add_argument('-vx', '--vx', dest='VXref', type=str, required=True,
help = 'input reference velocity x file path')
parser.add_argument('-vy', '--vy', dest='VYref', type=str, required=True,
help = 'input reference velocity y file path')
parser.add_argument('-ssm', '--ssm', dest='SSM', type=str, required=True,
help = 'input stable surface mask file path')
parser.add_argument('-o','--out', dest='output', type=str, required=True,
help = 'Variable to output')
parser.add_argument('-e','--epsg', dest='epsg', type=int, default=3413,
help = 'EPSG code')
return parser.parse_args()
def v_error_cal(vx_error, vy_error):
vx = np.random.normal(0, vx_error, 1000000)
vy = np.random.normal(0, vy_error, 1000000)
v = np.sqrt(vx2 + vy2)
return np.std(v)
if name == 'main':
#####Parse command line
inps = cmdLineParse()
xds = xarray.open_dataset(inps.nc)
VX = xds['vx'].data
VY = xds['vy'].data
V = xds['v'].data
V_error = xds['v_error'].data
if xds.attrs['scene_pair_type'] == 'radar':
VR = xds['vr'].data
VA = xds['va'].data
M11 = xds['M11'].data * xds['M11'].dr_to_vr_factor
M12 = xds['M12'].data * xds['M12'].dr_to_vr_factor
VXP = xds['vxp'].data
VYP = xds['vyp'].data
VP = xds['vp'].data
VP_error = xds['vp_error'].data
if np.logical_not(np.isnan(xds['vx'].stable_shift)):
VX += xds['vx'].stable_shift
VY += xds['vy'].stable_shift
if xds.attrs['scene_pair_type'] == 'radar':
VR += xds['vr'].stable_shift
VA += xds['va'].stable_shift
VXP += xds['vxp'].stable_shift
VYP += xds['vyp'].stable_shift
ds = gdal.Open(inps.VXref)
band = ds.GetRasterBand(1)
tran = ds.GetGeoTransform()
xoff = int(round((np.min(xds['vx'].x.data)-tran[0]-tran[1]/2)/tran[1]))
yoff = int(round((np.max(xds['vx'].y.data)-tran[3]-tran[5]/2)/tran[5]))
xcount = VX.shape[1]
ycount = VX.shape[0]
pdb.set_trace()
VXref = band.ReadAsArray(xoff, yoff, xcount, ycount)
ds = None
band = None
ds = gdal.Open(inps.VYref)
band = ds.GetRasterBand(1)
VYref = band.ReadAsArray(xoff, yoff, xcount, ycount)
ds = None
band = None
ds = gdal.Open(inps.SSM)
band = ds.GetRasterBand(1)
SSM = band.ReadAsArray(xoff, yoff, xcount, ycount)
SSM = SSM.astype('bool')
ds = None
band = None
pdb.set_trace()
if xds.attrs['scene_pair_type'] == 'radar':
VRref = M11 * VXref + M12 * VYref
VAref = VRref * 0.0
NoDataValue = -32767
VX and VY stable shift
stable_count = np.sum(SSM & np.logical_not(np.isnan(VX)))
V_temp = np.sqrt(VXref**2 + VYref**2)
try:
V_temp_threshold = np.percentile(V_temp[np.logical_not(np.isnan(V_temp))],25)
SSM1 = (V_temp <= V_temp_threshold)
except IndexError:
SSM1 = np.zeros(V_temp.shape).astype('bool')
stable_count1 = np.sum(SSM1 & np.logical_not(np.isnan(VX)))
vx_mean_shift = 0.0
vy_mean_shift = 0.0
vx_mean_shift1 = 0.0
vy_mean_shift1 = 0.0
if stable_count != 0:
temp = VX.copy() - VXref.copy()
temp[np.logical_not(SSM)] = np.nan
vx_mean_shift = np.median(temp[np.logical_not(np.isnan(temp))])-0.0
temp = VY.copy() - VYref.copy()
temp[np.logical_not(SSM)] = np.nan
vy_mean_shift = np.median(temp[np.logical_not(np.isnan(temp))])-0.0
if stable_count1 != 0:
temp = VX.copy() - VXref.copy()
temp[np.logical_not(SSM1)] = np.nan
vx_mean_shift1 = np.median(temp[np.logical_not(np.isnan(temp))])-0.0
temp = VY.copy() - VYref.copy()
temp[np.logical_not(SSM1)] = np.nan
vy_mean_shift1 = np.median(temp[np.logical_not(np.isnan(temp))])-0.0
if stable_count == 0:
if stable_count1 == 0:
stable_shift_applied = 0
else:
stable_shift_applied = 2
VX = VX - vx_mean_shift1
VY = VY - vy_mean_shift1
else:
stable_shift_applied = 1
VX = VX - vx_mean_shift
VY = VY - vy_mean_shift
VX and VY error
if stable_count != 0:
temp = VX.copy() - VXref.copy()
temp[np.logical_not(SSM)] = np.nan
vx_error_mask = np.std(temp[np.logical_not(np.isnan(temp))])
else:
vx_error_mask = np.nan
if stable_count1 != 0:
temp = VX.copy() - VXref.copy()
temp[np.logical_not(SSM1)] = np.nan
vx_error_slow = np.std(temp[np.logical_not(np.isnan(temp))])
else:
vx_error_slow = np.nan
if stable_shift_applied == 1:
vx_error = vx_error_mask
elif stable_shift_applied == 2:
vx_error = vx_error_slow
else:
vx_error = xds['vx'].vx_error_modeled
if stable_count != 0:
temp = VY.copy() - VYref.copy()
temp[np.logical_not(SSM)] = np.nan
vy_error_mask = np.std(temp[np.logical_not(np.isnan(temp))])
else:
vy_error_mask = np.nan
if stable_count1 != 0:
temp = VY.copy() - VYref.copy()
temp[np.logical_not(SSM1)] = np.nan
vy_error_slow = np.std(temp[np.logical_not(np.isnan(temp))])
else:
vy_error_slow = np.nan
if stable_shift_applied == 1:
vy_error = vy_error_mask
elif stable_shift_applied == 2:
vy_error = vy_error_slow
else:
vy_error = xds['vy'].vy_error_modeled
V and V_error
V = np.sqrt(VX**2+VY**2)
V_error = np.sqrt((vx_error * VX / V)**2 + (vy_error * VY / V)**2)
VX[np.isnan(VX)] = NoDataValue
VY[np.isnan(VY)] = NoDataValue
V[np.isnan(V)] = NoDataValue
VX = np.round(np.clip(VX, -32768, 32767)).astype(np.int16)
VY = np.round(np.clip(VY, -32768, 32767)).astype(np.int16)
V = np.round(np.clip(V, -32768, 32767)).astype(np.int16)
v_error = v_error_cal(vx_error, vy_error)
V_error[V==0] = v_error
V_error[np.isnan(V_error)] = NoDataValue
V_error = np.round(np.clip(V_error, -32768, 32767)).astype(np.int16)
Update nc file
xds['vx'].attrs['vx_error'] = int(round(vx_error*10))/10
if stable_shift_applied == 2:
xds['vx'].attrs['stable_shift'] = int(round(vx_mean_shift1*10))/10
elif stable_shift_applied == 1:
xds['vx'].attrs['stable_shift'] = int(round(vx_mean_shift*10))/10
else:
xds['vx'].attrs['stable_shift'] = np.nan
xds['vx'].attrs['stable_count_mask'] = stable_count
xds['vx'].attrs['stable_count_slow'] = stable_count1
if stable_count != 0:
xds['vx'].attrs['stable_shift_mask'] = int(round(vx_mean_shift10))/10
else:
xds['vx'].attrs['stable_shift_mask'] = np.nan
if stable_count1 != 0:
xds['vx'].attrs['stable_shift_slow'] = int(round(vx_mean_shift110))/10
else:
xds['vx'].attrs['stable_shift_slow'] = np.nan
if stable_count != 0:
xds['vx'].attrs['vx_error_mask'] = int(round(vx_error_mask*10))/10
else:
xds['vx'].attrs['vx_error_mask'] = np.nan
if stable_count1 != 0:
xds['vx'].attrs['vx_error_slow'] = int(round(vx_error_slow*10))/10
else:
xds['vx'].attrs['vx_error_slow'] = np.nan
xds['vy'].attrs['vy_error'] = int(round(vy_error10))/10
if stable_shift_applied == 2:
xds['vy'].attrs['stable_shift'] = int(round(vy_mean_shift110))/10
elif stable_shift_applied == 1:
xds['vy'].attrs['stable_shift'] = int(round(vy_mean_shift*10))/10
else:
xds['vy'].attrs['stable_shift'] = np.nan
xds['vy'].attrs['stable_count_mask'] = stable_count
xds['vy'].attrs['stable_count_slow'] = stable_count1
if stable_count != 0:
xds['vy'].attrs['stable_shift_mask'] = int(round(vy_mean_shift*10))/10
else:
xds['vy'].attrs['stable_shift_mask'] = np.nan
if stable_count1 != 0:
xds['vy'].attrs['stable_shift_slow'] = int(round(vy_mean_shift1*10))/10
else:
xds['vy'].attrs['stable_shift_slow'] = np.nan
if stable_count != 0:
xds['vy'].attrs['vy_error_mask'] = int(round(vy_error_mask10))/10
else:
xds['vy'].attrs['vy_error_mask'] = np.nan
if stable_count1 != 0:
xds['vy'].attrs['vy_error_slow'] = int(round(vy_error_slow10))/10
else:
xds['vy'].attrs['vy_error_slow'] = np.nan
xds['vx'].data = VX
xds['vy'].data = VY
xds['v'].data = V
xds['v_error'].data = V_error
if xds.attrs['scene_pair_type'] == 'radar':
# VR and VA stable shift
vr_mean_shift = 0.0
va_mean_shift = 0.0
vr_mean_shift1 = 0.0
va_mean_shift1 = 0.0
if stable_count != 0:
temp = VR.copy() - VRref.copy()
temp[np.logical_not(SSM)] = np.nan
vr_mean_shift = np.median(temp[np.logical_not(np.isnan(temp))])-0.0
temp = VA.copy() - VAref.copy()
temp[np.logical_not(SSM)] = np.nan
va_mean_shift = np.median(temp[np.logical_not(np.isnan(temp))])-0.0
if stable_count1 != 0:
temp = VR.copy() - VRref.copy()
temp[np.logical_not(SSM1)] = np.nan
vr_mean_shift1 = np.median(temp[np.logical_not(np.isnan(temp))])-0.0
temp = VA.copy() - VAref.copy()
temp[np.logical_not(SSM1)] = np.nan
va_mean_shift1 = np.median(temp[np.logical_not(np.isnan(temp))])-0.0
if stable_count == 0:
if stable_count1 == 0:
stable_shift_applied = 0
else:
stable_shift_applied = 2
VR = VR - vr_mean_shift1
VA = VA - va_mean_shift1
else:
stable_shift_applied = 1
VR = VR - vr_mean_shift
VA = VA - va_mean_shift
# VR and VA error
if stable_count != 0:
temp = VR.copy() - VRref.copy()
temp[np.logical_not(SSM)] = np.nan
vr_error_mask = np.std(temp[np.logical_not(np.isnan(temp))])
else:
vr_error_mask = np.nan
if stable_count1 != 0:
temp = VR.copy() - VRref.copy()
temp[np.logical_not(SSM1)] = np.nan
vr_error_slow = np.std(temp[np.logical_not(np.isnan(temp))])
else:
vr_error_slow = np.nan
if stable_shift_applied == 1:
vr_error = vr_error_mask
elif stable_shift_applied == 2:
vr_error = vr_error_slow
else:
vr_error = xds['vr'].vr_error_modeled
if stable_count != 0:
temp = VA.copy() - VAref.copy()
temp[np.logical_not(SSM)] = np.nan
va_error_mask = np.std(temp[np.logical_not(np.isnan(temp))])
else:
va_error_mask = np.nan
if stable_count1 != 0:
temp = VA.copy() - VAref.copy()
temp[np.logical_not(SSM1)] = np.nan
va_error_slow = np.std(temp[np.logical_not(np.isnan(temp))])
else:
va_error_slow = np.nan
if stable_shift_applied == 1:
va_error = va_error_mask
elif stable_shift_applied == 2:
va_error = va_error_slow
else:
va_error = xds['va'].va_error_modeled
# V and V_error
VR[np.isnan(VXP)] = NoDataValue
VA[np.isnan(VYP)] = NoDataValue
VR = np.round(np.clip(VR, -32768, 32767)).astype(np.int16)
VA = np.round(np.clip(VA, -32768, 32767)).astype(np.int16)
# Update nc file
xds['vr'].attrs['vr_error'] = int(round(vr_error*10))/10
if stable_shift_applied == 2: