203
203
dmin = args .radius [0 ] * 1e3 # min search radius (km -> m)
204
204
dmax = args .radius [1 ] * 1e3 + 1e-4 # max search radius (km -> m)
205
205
dres_ = args .resparam [0 ] # resolution param for weighting func [1]
206
- nreloc = args .nreloc [0 ] # number of relocations
206
+ nreloc = args .nreloc [0 ] # number of relocations
207
207
nlim = args .minobs [0 ] # min obs for solution
208
208
niter = args .niter [0 ] # number of iterations for solution
209
209
tspan = args .tspan # min/max time for solution (d.yr)
229
229
# [4] If err and id cols = -1, then they are not used.
230
230
231
231
232
- def binning (x , y , xmin = None , xmax = None , dx = 1 / 12. ,
232
+ def binning (x , y , xmin = None , xmax = None , dx = 1 / 12. ,
233
233
window = 3 / 12. , interp = False , median = False ):
234
234
"""Time-series binning (w/overlapping windows).
235
235
@@ -304,7 +304,7 @@ def sigma_filter(x, y, order=1, window=3/12., n_iter=3, n_sigma=3):
304
304
if len (idx ) == 0 : break # if no data to filter, stop iterating
305
305
y_res [idx ] = np .nan
306
306
if np .sum (~ np .isnan (y_res )) < 10 : break ##NOTE: Arbitrary min obs
307
- y_filt [np .isnan (y_res )] = np .nan
307
+ y_filt [np .isnan (y_res )] = np .nan
308
308
return y_filt
309
309
310
310
@@ -360,7 +360,7 @@ def get_radius_idx(x, y, x0, y0, r, Tree, n_reloc=0):
360
360
# Either no relocation or not enough points to do relocation
361
361
if n_reloc < 1 or len (idx ) < 2 : return idx , reloc_dist
362
362
363
- # Relocate center of search radius and query again
363
+ # Relocate center of search radius and query again
364
364
for k in range (n_reloc ):
365
365
366
366
# Compute new search location => relocate initial center
@@ -400,7 +400,7 @@ def is_empty(ifile):
400
400
401
401
# Main function for computing parameters
402
402
def main (ifile , n = '' , robust_fit = True , n_iter = niter ):
403
-
403
+
404
404
# Check for empty file
405
405
if is_empty (ifile ):
406
406
print (('SKIP FILE: EMPTY OR CORRUPTED FILE:' , ifile ))
@@ -471,9 +471,9 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
471
471
472
472
# Check bbox for obs.
473
473
if len (x [Ig ]) == 0 :
474
- print (('SKIP FILE: NO DATA POINTS INSIDE BBOX:' , ifile ))
474
+ print (('SKIP FILE: NO DATA POINTS INSIDE BBOX:' , ifile ))
475
475
return
476
-
476
+
477
477
print (('Number of obs. edited by bbox!' , 'before:' , len (x ), 'after:' , len (x [Ig ])))
478
478
479
479
# Only select wanted data
@@ -519,7 +519,7 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
519
519
# Grid solution - defined by nodes
520
520
Xi , Yi = make_grid (xmin , xmax , ymin , ymax , dx , dy )
521
521
522
- xi , yi = Xi .ravel (), Yi .ravel ()
522
+ xi , yi = Xi .ravel (), Yi .ravel ()
523
523
coord = list (zip (x .ravel (), y .ravel ()))
524
524
525
525
print ('building the k-d tree ...' )
@@ -563,11 +563,11 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
563
563
564
564
if len (i_cell ) < nlim : continue # use larger radius
565
565
566
- tcap , hcap = time [i_cell ], height [i_cell ]
566
+ tcap , hcap = time [i_cell ], height [i_cell ]
567
567
568
568
Nb = sum (~ np .isnan (hcap )) # length before editing
569
569
570
- # 3-sigma filter
570
+ # 3-sigma filter
571
571
if SIGMAFILT :
572
572
#hcap = sigma_filter(tcap, hcap, order=1, n_sigma=3, n_iter=3) ##NOTE: It removes too much!!!
573
573
hcap [np .abs ( hcap - np .nanmedian (hcap ) ) > mad_std (hcap ) * 3 ] = np .nan
@@ -578,7 +578,7 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
578
578
n_mon , t_span = n_months (tcap , hcap , tstep = tstep )
579
579
580
580
##NOTE: Not using n_mon and t_span to constrain the solution! <<<<<<<<<<<<<<<<<<<<<
581
- # If enough data accept radius
581
+ # If enough data accept radius
582
582
#if Na >= nlim and n_mon >= MINMONTHS and t_span >= dtlim:
583
583
if Na >= nlim :
584
584
break
@@ -609,8 +609,8 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
609
609
xc = np .median (xcap ) # update inversion cell coords
610
610
yc = np .median (ycap )
611
611
612
- # Define resolution param (a fraction of the accepted radius)
613
- dres = dres_ * rad
612
+ # Define resolution param (a fraction of the accepted radius)
613
+ dres = dres_ * rad
614
614
615
615
# Estimate variance
616
616
vcap = scap * scap
@@ -622,7 +622,7 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
622
622
tref = np .nanmean (tcap )
623
623
else :
624
624
tref = np .float (tref_ )
625
-
625
+
626
626
# Design matrix elements
627
627
c0 = np .ones (len (xcap )) # intercept (0)
628
628
c1 = xcap - xc # dx (1)
@@ -645,7 +645,7 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
645
645
Wcap = 1.0 / (vcap * (1.0 + (dist / dres )* (dist / dres )))
646
646
647
647
# Create some intermediate output variables
648
- sx , sy , at , ae , bi = np .nan , np .nan , np .nan , np .nan , np .nan
648
+ sx , sy , at , ae , bi = np .nan , np .nan , np .nan , np .nan , np .nan
649
649
650
650
# Setup design matrix
651
651
if model == 0 :
@@ -666,7 +666,7 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
666
666
mcol = [6 , 7 , 8 , 9 ]
667
667
668
668
has_bias = False # bias flag
669
-
669
+
670
670
# Check if bias is needed
671
671
if len (np .unique (mcap )) > 1 :
672
672
# Add bias to design matrix
@@ -694,10 +694,10 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
694
694
Cm = model_fit .params # coeffs
695
695
Ce = model_fit .bse # std err
696
696
resid = model_fit .resid # data - model
697
-
697
+
698
698
# Check rate and error
699
699
if np .abs (Cm [- 1 ]) > dhlim or np .isinf (Ce [- 1 ]): continue ##NOTE: Important for ICESat !!!
700
-
700
+
701
701
# Residuals dH = H - A * Cm (remove linear trend)
702
702
dh = hcap - np .dot (Acap , Cm )
703
703
@@ -796,7 +796,7 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
796
796
DATA2 = np .delete (DATA2 .T , i_nan , 1 ).T
797
797
else :
798
798
##NOTE: NaNs are not removed in case a grid soluction (n_reloc=0) is selected.
799
- if not nreloc : grids = [d .reshape (Xi .shape ) for d in DATA0 .T ] # 1d -> 2d (grids)
799
+ if not nreloc : grids = [d .reshape (Xi .shape ) for d in DATA0 .T ] # 1d -> 2d (grids)
800
800
801
801
variables = ['lat' , 'lon' , 'trend' , 'trend_err' , 'accel' , 'accel_err' ,
802
802
'height' , 'height_err' , 'model_rms' , 'slope_x' , 'slope_y' ,
@@ -830,7 +830,7 @@ def main(ifile, n='', robust_fit=True, n_iter=niter):
830
830
for v ,a in zip (variables , DATA0 .T ): fo0 [v ] = a # 1d arrays
831
831
else :
832
832
for v ,g in zip (variables , grids ): fo0 [v ] = g # 2d arrays
833
- fo0 ['x' ], fo0 ['y' ] = Xi [0 ,:], Yi [:,0 ]
833
+ fo0 ['x' ], fo0 ['y' ] = Xi [0 ,:], Yi [:,0 ]
834
834
835
835
# Save binned time series values
836
836
with h5py .File (ofile1 , 'w' ) as fo1 :
0 commit comments