forked from insarlab/MintPy
-
Notifications
You must be signed in to change notification settings - Fork 0
/
utils1.py
952 lines (818 loc) · 36.9 KB
/
utils1.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
"""Miscellaneous utilities - dependent on utils0."""
############################################################
# Program is part of MintPy #
# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi #
# Author: Zhang Yunjun, Heresh Fattahi, 2013 #
############################################################
# Recommend import:
# from mintpy.utils import utils as ut
import glob
import os
import re
import shutil
import time
import h5py
import numpy as np
import mintpy
from mintpy.objects import GEOMETRY_DSET_NAMES, deramp, ifgramStack, timeseries
from mintpy.utils import ptime, readfile, writefile
from mintpy.utils.utils0 import *
#################################### Geometry #########################################
def get_center_lat_lon(geom_file, box=None):
"""Get the lat/lon of the scene center"""
meta = readfile.read_attribute(geom_file)
if box is None:
box = (0, 0, int(meta['WIDTH']), int(meta['LENGTH']))
col_c = int((box[0] + box[2]) / 2)
row_c = int((box[1] + box[3]) / 2)
if 'Y_FIRST' in meta.keys():
lat0 = float(meta['Y_FIRST'])
lon0 = float(meta['X_FIRST'])
lat_step = float(meta['Y_STEP'])
lon_step = float(meta['X_STEP'])
lat_c = lat0 + lat_step * row_c
lon_c = lon0 + lon_step * col_c
else:
box_c = (col_c, row_c, col_c+1, row_c+1)
lat_c = float(readfile.read(geom_file, datasetName='latitude', box=box_c)[0])
lon_c = float(readfile.read(geom_file, datasetName='longitude', box=box_c)[0])
return lat_c, lon_c
#################################### Data Operation ###################################
def get_residual_std(timeseries_resid_file, mask_file='maskTempCoh.h5', ramp_type='quadratic'):
"""Calculate deramped standard deviation in space for each epoch of input timeseries file.
Parameters: timeseries_resid_file - string, timeseries HDF5 file,
e.g. timeseries_ERA5_demErrInvResid.h5
mask_file - string, mask file, e.g. maskTempCoh.h5
ramp_type - string, ramp type, e.g. linear, quadratic, no for do not remove ramp
Returns: std_list - list of float, standard deviation of deramped input timeseries file
date_list - list of string in YYYYMMDD format, corresponding dates
std_file - string, text file with std and date info.
Example: import mintpy.utils.utils as ut
std_list, date_list = ut.get_residual_std('timeseries_ERA5_demErrInvResid.h5',
'maskTempCoh.h5')[:2]
"""
# Intermediate files name
# ramp_type can sometimes be False, thus, should be treated the same as "no"
if not ramp_type or ramp_type == 'no':
print('No ramp removal')
deramped_file = timeseries_resid_file
else:
deramped_file = f'{os.path.splitext(timeseries_resid_file)[0]}_ramp.h5'
std_file = os.path.splitext(deramped_file)[0]+'_std.txt'
# Get residual std text file
if run_or_skip(out_file=std_file, in_file=[timeseries_resid_file, mask_file], readable=False) == 'run':
if run_or_skip(out_file=deramped_file, in_file=timeseries_resid_file) == 'run':
if not os.path.isfile(timeseries_resid_file):
msg = 'Can not find input timeseries residual file: '+timeseries_resid_file
msg += '\nRe-run dem_error.py to generate it.'
raise Exception(msg)
else:
#print('removing a {} ramp from file: {}'.format(ramp_type, timeseries_resid_file))
deramped_file = run_deramp(
timeseries_resid_file,
ramp_type=ramp_type,
mask_file=mask_file,
out_file=deramped_file,
)
print('calculating residual standard deviation for each epoch from file: '+deramped_file)
std_file = timeseries(deramped_file).timeseries_std(maskFile=mask_file, outFile=std_file)
# Read residual std text file
print('read timeseries RMS from file: '+std_file)
fc = np.loadtxt(std_file, dtype=bytes).astype(str)
std_list = fc[:, 1].astype(np.float32).tolist()
date_list = list(fc[:, 0])
return std_list, date_list, std_file
def get_residual_rms(timeseries_resid_file, mask_file='maskTempCoh.h5', ramp_type='quadratic'):
"""Calculate deramped Root Mean Square in space for each epoch of input timeseries file.
Parameters: timeseries_resid_file : string,
timeseries HDF5 file, e.g. timeseries_ERA5_demErrInvResid.h5
mask_file : string,
mask file, e.g. maskTempCoh.h5
ramp_type : string,
ramp type, e.g. linear, quadratic, no for do not remove ramp
Returns: rms_list : list of float,
Root Mean Square of deramped input timeseries file
date_list : list of string in YYYYMMDD format,
corresponding dates
rms_file : string, text file with rms and date info.
Example:
import mintpy.utils.utils as ut
rms_list, date_list = ut.get_residual_rms('timeseriesResidual.h5', 'maskTempCoh.h5')
"""
# Intermediate files name
# ramp_type can sometimes be False, thus, should be treated the same as "no"
if not ramp_type or ramp_type == 'no':
print('No ramp removal')
deramped_file = timeseries_resid_file
else:
deramped_file = f'{os.path.splitext(timeseries_resid_file)[0]}_ramp.h5'
fdir = os.path.dirname(os.path.abspath(deramped_file))
fbase = os.path.splitext(os.path.basename(deramped_file))[0]
rms_file = os.path.join(fdir, f'rms_{fbase}.txt')
# Get residual RMS text file
if run_or_skip(out_file=rms_file, in_file=[timeseries_resid_file, mask_file], readable=False) == 'run':
if run_or_skip(out_file=deramped_file, in_file=timeseries_resid_file) == 'run':
if not os.path.isfile(timeseries_resid_file):
msg = 'Can not find input timeseries residual file: '+timeseries_resid_file
msg += '\nRe-run dem_error.py to generate it.'
raise Exception(msg)
else:
#print('remove {} ramp from file: {}'.format(ramp_type, timeseries_resid_file))
deramped_file = run_deramp(
timeseries_resid_file,
ramp_type=ramp_type,
mask_file=mask_file,
out_file=deramped_file,
)
print('\ncalculating residual RMS for each epoch from file: '+deramped_file)
rms_file = timeseries(deramped_file).timeseries_rms(
maskFile=mask_file,
outFile=rms_file,
)
# Read residual RMS text file
print('read timeseries residual RMS from file: '+rms_file)
fc = np.loadtxt(rms_file, dtype=bytes).astype(str)
rms_list = fc[:, 1].astype(np.float32).tolist()
date_list = list(fc[:, 0])
return rms_list, date_list, rms_file
def nonzero_mask(File, out_file='maskConnComp.h5', datasetName=None):
"""Generate mask file for non-zero value of input multi-group hdf5 file"""
atr = readfile.read_attribute(File)
k = atr['FILE_TYPE']
if k == 'ifgramStack':
mask = ifgramStack(File).nonzero_mask(datasetName=datasetName)
else:
print('Only ifgramStack file is supported for now, input is '+k)
return None
atr['FILE_TYPE'] = 'mask'
writefile.write(mask, out_file=out_file, metadata=atr)
return out_file
def spatial_average(fname, datasetName='coherence', maskFile=None, box=None,
saveList=False, checkAoi=True, reverseMask=False, threshold=None):
"""Read/Calculate Spatial Average of input file.
If input file is text file, read it directly;
If input file is data matrix file:
If corresponding text file exists with the same mask file/AOI info, read it directly;
Otherwise, calculate it from data file.
Only non-nan pixel is considered.
Parameters: fname - string, path of input file
maskFile - string, path of mask file, e.g. maskTempCoh.h5
box - 4-tuple defining the left, upper, right, and lower pixel coordinate
saveList - bool, save (list of) mean value into text file
reverseMask - bool, perform analysis within masked regions instead of outside of them
threshold - float, calculate area ratio above threshold instead of spatial average
Returns: meanList - list(float) or float, average value in space for each epoch of input file
dateList - list(str) or str, for date info
date12_list, e.g. 101120-110220, for interferograms/coherence
date8_list, e.g. 20101120, for timeseries
file name, e.g. velocity.h5, for all the other file types
Example: meanList = spatial_average('inputs/ifgramStack.h5')[0]
meanList, date12_list = spatial_average('inputs/ifgramStack.h5',
maskFile='maskTempCoh.h5',
saveList=True)
"""
def read_text_file(fname):
txtContent = np.loadtxt(fname, dtype=bytes).astype(str)
meanList = [float(i) for i in txtContent[:, 1]]
dateList = [i for i in txtContent[:, 0]]
return meanList, dateList
# Baic File Info
atr = readfile.read_attribute(fname)
k = atr['FILE_TYPE']
if not box:
box = (0, 0, int(atr['WIDTH']), int(atr['LENGTH']))
# default output filename
prefix = datasetName if k == 'ifgramStack' else os.path.splitext(os.path.basename(fname))[0]
suffix = 'SpatialAvg' if threshold is None else 'AreaRatio'
suffix += 'RevMsk' if reverseMask else ''
txtFile = prefix + suffix + '.txt'
# If input is text file
if fname.endswith(suffix):
print('Input file is spatial average txt already, read it directly')
meanList, dateList = read_text_file(fname)
return meanList, dateList
# Read existing txt file only if 1) data file is older AND 2) same AOI
file_line = f'# Data file: {os.path.basename(fname)}\n'
mask_line = f'# Mask file: {maskFile}\n'
aoi_line = f'# AOI box: {box}\n'
thres_line = f'# Threshold: {threshold}\n'
try:
# Read AOI line from existing txt file
fl = open(txtFile)
lines = fl.readlines()
fl.close()
# 1. aoi
if checkAoi:
try:
aoi_line_orig = [i for i in lines if '# AOI box:' in i][0]
except:
aoi_line_orig = ''
else:
aoi_line_orig = aoi_line
# 2. mask file
try:
mask_line_orig = [i for i in lines if '# Mask file:' in i][0]
except:
mask_line_orig = ''
# 3. mask file - modification time
update_mask_file = run_or_skip(out_file=txtFile, in_file=[maskFile], readable=False)
# 4. data file - modification time
if k == 'ifgramStack':
with h5py.File(fname, 'r') as f:
ti = float(f[datasetName].attrs.get('MODIFICATION_TIME', os.path.getmtime(fname)))
else:
ti = os.path.getmtime(fname)
to = os.path.getmtime(txtFile)
if (aoi_line_orig == aoi_line
and mask_line_orig == mask_line
and update_mask_file == 'skip'
and ti <= to):
print(txtFile+' already exists, read it directly')
meanList, dateList = read_text_file(txtFile)
return meanList, dateList
except:
pass
# use median instead of mean for offset measurement
if datasetName and 'offset' in datasetName:
useMedian = True
else:
useMedian = False
# Calculate mean coherence or area ratio list
if k == 'ifgramStack':
obj = ifgramStack(fname)
obj.open(print_msg=False)
meanList, dateList = obj.spatial_average(
datasetName=datasetName,
maskFile=maskFile,
box=box,
useMedian=useMedian,
reverseMask=reverseMask,
threshold=threshold,
)
pbase = obj.pbaseIfgram
tbase = obj.tbaseIfgram
obj.close()
elif k == 'timeseries':
meanList, dateList = timeseries(fname).spatial_average(
maskFile=maskFile,
box=box,
reverseMask=reverseMask,
threshold=threshold,
)
else:
data = readfile.read(fname, box=box)[0]
if maskFile and os.path.isfile(maskFile):
print('mask from file: '+maskFile)
mask = readfile.read(maskFile, datasetName='mask', box=box)[0]
data[mask == int(reverseMask)] = np.nan
# calculate area ratio if threshold is specified
# percentage of pixels with value above the threshold
if threshold is not None:
data[data > threshold] = 1
data[data <= threshold] = 0
meanList = np.nanmean(data)
dateList = [os.path.basename(fname)]
# Write mean coherence list into text file
if saveList:
print('write average value in space into text file: '+txtFile)
fl = open(txtFile, 'w')
# Write comments
fl.write(file_line+mask_line+aoi_line+thres_line)
# Write data list
numLine = len(dateList)
if k == 'ifgramStack':
fl.write('#\tDATE12\t\tMean\tBtemp/days\tBperp/m\t\tNum\n')
for i in range(numLine):
fl.write('%s\t%.4f\t%8.0f\t%8.1f\t%d\n' %
(dateList[i], meanList[i], tbase[i], pbase[i], i))
else:
fl.write('#\tDATE12\t\tMean\n')
for i in range(numLine):
fl.write(f'{dateList[i]}\t{meanList[i]:.4f}\n')
fl.close()
# read from text file (in 1e-4 precision)
# to ensure output value consistency
meanList, dateList = read_text_file(txtFile)
if len(meanList) == 1:
meanList = meanList[0]
dateList = dateList[0]
return meanList, dateList
def temporal_average(fname, datasetName='coherence', updateMode=False, outFile=None):
"""Calculate temporal average of multi-temporal dataset, equivalent to stacking
For ifgramStack/unwrapPhase, return average phase velocity
Parameters: fname - str, file to be averaged in time
datasetName - str, dataset to be read from input file, for multiple
datasets file - ifgramStack - only
e.g.: coherence, unwrapPhase
updateMode - bool
outFile - str, output filename
None for auto output filename
False for do not save as output file
Returns: dataMean - 2D np.ndarray
outFile - str, output file name
Examples: avgPhaseVel = ut.temporal_average('ifgramStack.h5', datasetName='unwrapPhase')[0]
ut.temporal_average('ifgramStack.h5', datasetName='coherence',
outFile='avgSpatialCoh.h5', updateMode=True)
"""
atr = readfile.read_attribute(fname, datasetName=datasetName)
k = atr['FILE_TYPE']
if k not in ['ifgramStack', 'timeseries']:
print(f'WARNING: input file is not multi-temporal file: {fname}, return itself.')
data = readfile.read(fname)[0]
return data, fname
# Default output filename
if outFile is None:
ext = os.path.splitext(fname)[1]
if not outFile:
if k == 'ifgramStack':
if datasetName == 'coherence':
outFile = 'avgSpatialCoh.h5'
elif 'unwrapPhase' in datasetName:
outFile = 'avgPhaseVelocity.h5'
else:
outFile = f'avg{datasetName}.h5'
elif k == 'timeseries':
if k in fname:
processMark = os.path.basename(fname).split('timeseries')[1].split(ext)[0]
outFile = f'avgDisplacement{processMark}.h5'
else:
outFile = f'avg{fname}.h5'
if updateMode and os.path.isfile(outFile):
dataMean = readfile.read(outFile)[0]
return dataMean, outFile
# Calculate temporal average
if k == 'ifgramStack':
dataMean = ifgramStack(fname).temporal_average(datasetName=datasetName)
if 'unwrapPhase' in datasetName:
atr['FILE_TYPE'] = 'velocity'
atr['UNIT'] = 'm/year'
else:
atr['FILE_TYPE'] = datasetName
elif k == 'timeseries':
dataMean = timeseries(fname).temporal_average()
atr['FILE_TYPE'] = 'displacement'
if outFile:
writefile.write(dataMean, out_file=outFile, metadata=atr)
return dataMean, outFile
#################################### File IO ##########################################
def get_file_list(file_list, abspath=False, coord=None):
"""Get all existed files matching the input list of file pattern
Parameters: file_list - string or list of string, input file/directory pattern
abspath - bool, return absolute path or not
coord - string, return files with specific coordinate type: geo or radar
if none, skip the checking and return all files
Returns: file_list_out - list of string, existed file path/name, [] if not existed
Example: file_list = get_file_list(['*velocity*.h5','timeseries*.h5'])
file_list = get_file_list('timeseries*.h5')
"""
if not file_list:
return []
if isinstance(file_list, str):
file_list = [file_list]
# Get rid of None element
file_list = [x for x in file_list if x is not None]
file_list_out = []
for fname in file_list:
fnames = glob.glob(fname)
file_list_out += sorted(list(set(fnames) - set(file_list_out)))
if abspath:
file_list_out = [os.path.abspath(i) for i in file_list_out]
if coord is not None:
for fname in list(file_list_out):
atr = readfile.read_attribute(fname)
if coord in ['geo']:
if 'Y_FIRST' not in atr.keys():
file_list_out.remove(fname)
elif coord in ['radar', 'rdr', 'rdc']:
if 'Y_FIRST' in atr.keys():
file_list_out.remove(fname)
else:
msg = f'un-recognized input coord type: {coord}'
raise ValueError(msg)
return file_list_out
def get_lookup_file(filePattern=None, abspath=False, print_msg=True):
"""Find lookup table file with/without input file pattern
Parameters: filePattern - list of str
abspath - bool, return absolute path or not
print_msg - bool, printout message or not
Returns: outFile - str, path of the lookup file
"""
# Search Existing Files
if not filePattern:
fileList = ['geometryRadar.h5',
'geometryGeo_tight.h5', 'geometryGeo.h5',
'geomap*lks_tight.trans', 'geomap*lks.trans',
'sim*_tight.UTM_TO_RDC', 'sim*.UTM_TO_RDC']
dirList = ['inputs', '', '../inputs']
# file/dirList --> filePattern
filePattern = []
for dirname in dirList:
filePattern += [os.path.join(dirname, fname) for fname in fileList]
existFiles = []
try:
existFiles = get_file_list(filePattern)
except:
if print_msg:
print('ERROR: No geometry / lookup table file found!')
print('It should be like:')
print(filePattern)
return None
# Check Files Info
outFile = None
for fname in existFiles:
readfile.read_attribute(fname)
for dsName in ['longitude', 'rangeCoord']:
try:
readfile.read(fname, datasetName=dsName, print_msg=False)
outFile = fname
break
except:
pass
if not outFile:
if print_msg:
print('No lookup table (longitude or rangeCoord) found in files.')
return None
# Path Format
if abspath:
outFile = os.path.abspath(outFile)
return outFile
def get_geometry_file(dset_list, work_dir=None, coord='geo', abspath=True, print_msg=True):
"""Find geometry file containing input specific dataset"""
if isinstance(dset_list, str):
dset_list = [dset_list]
for dset in dset_list:
if dset not in GEOMETRY_DSET_NAMES:
raise ValueError(f'unrecognized geometry dataset name: {dset}')
if not work_dir:
work_dir = os.getcwd()
# search *geometry*.h5 files
fname_list = [os.path.join(work_dir, i) for i in ['*geometry*.h5', '*/*geometry*.h5', '../*/geometry*.h5']]
fname_list = get_file_list(fname_list, coord=coord)
if len(fname_list) == 0:
if print_msg:
print('No geometry file found.')
return None
# check dset in the existing h5 files
for fname in list(fname_list): #use list() as temp copy to handle modifying list during the loop
if any(dset not in readfile.get_dataset_list(fname) for dset in dset_list):
fname_list.remove(fname)
if len(fname_list) == 0:
if print_msg:
print(f'No geometry file with dataset {dset_list} found')
return None
geom_file = fname_list[0]
if abspath:
geom_file = os.path.abspath(geom_file)
return geom_file
def update_template_file(template_file, extra_dict, delimiter='='):
"""Update option value in template_file with value from input extra_dict"""
# Compare and skip updating template_file if no new option value found.
update = False
orig_dict = readfile.read_template(template_file)
for key, value in orig_dict.items():
if key in extra_dict.keys() and extra_dict[key] != value:
update = True
if not update:
print('No new option value found, skip updating '+template_file)
return template_file
# Update template_file with new value from extra_dict
tmp_file = template_file+'.tmp'
f_tmp = open(tmp_file, 'w')
for line in open(template_file):
c = [i.strip() for i in line.strip().split(delimiter, 1)]
if not line.startswith(('%', '#')) and len(c) > 1:
key = c[0]
value = str.replace(c[1], '\n', '').split("#")[0].strip()
if key in extra_dict.keys() and extra_dict[key] != value:
# prepare value string to search & replace following "re" expression syntax
# link: https://docs.python.org/3/library/re.html
value2search = value
# 1. interpret special symbols as characters
for symbol in ['*', '[', ']', '(', ')']:
value2search = value2search.replace(symbol, fr"\{symbol}")
# 2. use "= {OLD_VALUE}" for search/replace to be more robust
# against the scenario when key name contains {OLD_VALUE}
# i.e. mintpy.load.autoPath
value2search = delimiter+r'[\s]*'+value2search
old_value_str = re.findall(value2search, line)[0]
new_value_str = old_value_str.replace(value, extra_dict[key])
line = line.replace(old_value_str, new_value_str, 1)
print(f' {key}: {value} --> {extra_dict[key]}')
f_tmp.write(line)
f_tmp.close()
# Overwrite existing original template file
shutil.move(tmp_file, template_file)
return template_file
def add_attribute(fname, atr_new=dict(), print_msg=False):
"""Add/update input attribute of the give file.
Parameters: fname - string, path/name of file
atr_new - dict, attributes to be added/updated
if value is None, delete the item from input file attributes
Returns: fname - string, path/name of updated file
"""
vprint = print if print_msg else lambda *args, **kwargs: None
# read existing attributes
atr = readfile.read_attribute(fname)
key_list = list(atr.keys())
# compare new attributes with existing ones
update = update_attribute_or_not(atr_new, atr)
if not update:
vprint('All updated (removed) attributes already exists (do not exists)'
' and have the same value, skip update.')
return fname
# update attributes in the input data file
fext = os.path.splitext(fname)[1]
if fext in ['.h5', '.he5']:
with h5py.File(fname, 'r+') as f:
for key, value in iter(atr_new.items()):
if value == 'None' or value is None:
# delete the item for invalid input (None)
if key in key_list:
f.attrs.pop(key)
vprint(f'remove {key}')
else:
# update the item for valid input
f.attrs[key] = str(value)
vprint(f'add/update {key} = {str(value)}')
else:
for key, value in iter(atr_new.items()):
if value == 'None' or value is None:
# delete the item for invalid input (None)
if key in key_list:
atr.pop(key)
vprint(f'remove {key}')
else:
# update the item for valid input
atr[key] = str(value)
vprint(f'add/update {key} = {str(value)}')
# write to RSC file
writefile.write_roipac_rsc(atr, fname+'.rsc', print_msg=print_msg)
return fname
def check_file_size(fname_list, mode_width=None, mode_length=None):
"""Check file size in the list of files, and drop those not in the same size with majority."""
# If input file list is empty
if not fname_list:
return fname_list, None, None
# Read Width/Length list
width_list = []
length_list = []
for fname in fname_list:
atr = readfile.read_attribute(fname)
width_list.append(atr['WIDTH'])
length_list.append(atr['LENGTH'])
# Mode of Width and Length
mode_width = mode_width if mode_width else most_common(width_list)
mode_length = mode_length if mode_length else most_common(length_list)
# Update Input List
fname_list_out = list(fname_list)
if (width_list.count(mode_width) != len(width_list)
or length_list.count(mode_length) != len(length_list)):
print('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
print('WARNING: Some files may have the wrong dimensions!')
print('All files should have the same size.')
print('The width and length of the majority of files are: %s, %s' %
(mode_width, mode_length))
print('But the following files have different dimensions and thus will not be loaded:')
for fname, length, width in zip(fname_list, length_list, width_list):
if width != mode_width or length != mode_length:
print(f'{fname} width: {width} length: {length}')
fname_list_out.remove(fname)
print('\nNumber of files left: '+str(len(fname_list_out)))
print('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
return fname_list_out, mode_width, mode_length
#################################### Interaction ##########################################
def is_file_exist(file_list, abspath=True):
"""Check if any file in the file list 1) exists and 2) readable
Parameters: file_list : str or list(str), file name with/without wildcards
abspath : bool, return absolute file name/path or not
Returns: file_path : string, found file name/path; None if not.
"""
try:
file = get_file_list(file_list, abspath=abspath)[0]
readfile.read_attribute(file)
except:
file = None
return file
def run_or_skip(out_file, in_file=None, readable=True, print_msg=True):
"""Check whether to update out_file or not.
return run if any of the following meets:
1. out_file is empty, e.g. None, []
2. out_file is not existed
3. out_file is not readable by readfile.read_attribute() when readable=True
4. out_file is older than in_file, if in_file is not None
Otherwise, return skip.
If in_file=None and out_file exists and readable, return skip
Parameters: out_file - string or list of string, output file(s)
in_file - string or list of string, input file(s)
readable - bool, check if the 1st output file has attribute 'WIDTH'
print_msg - bool, print message
Returns: run/skip - str, whether to update output file or not
Example: if ut.run_or_skip(out_file='timeseries_ERA5_demErr.h5', in_file='timeseries_ERA5.h5'):
if ut.run_or_skip(out_file='exclude_date.txt',
in_file=['timeseries_ERA5_demErrInvResid.h5',
'maskTempCoh.h5',
'smallbaselineApp.cfg'],
readable=False):
"""
# 1 - check existence of output files
if not out_file:
return 'run'
else:
if isinstance(out_file, str):
out_file = [out_file]
if not all(os.path.isfile(i) for i in out_file):
return 'run'
# 2 - check readability of output files
if readable:
try:
readfile.read_attribute(out_file[0])['WIDTH']
except:
if print_msg:
print(f'{out_file[0]} exists, but can not read, remove it.')
os.remove(out_file[0])
return 'run'
# 3 - check modification time of output and input files
if in_file:
in_file = get_file_list(in_file)
# Check modification time
if in_file:
t_in = max(os.path.getmtime(i) for i in in_file)
t_out = min(os.path.getmtime(i) for i in out_file)
if t_in > t_out:
return 'run'
elif print_msg:
print(f'{out_file} exists and is newer than {in_file} --> skip.')
return 'skip'
def check_template_auto_value(templateDict, auto_file='defaults/smallbaselineApp_auto.cfg'):
"""Replace auto value based on the input auto config file."""
## Read default template value and turn yes/no to True/False
templateAutoFile = os.path.join(os.path.dirname(mintpy.__file__), auto_file)
templateAutoDict = readfile.read_template(templateAutoFile)
# if cluster != local, change auto value of numWorker
cluster_key = 'mintpy.compute.cluster'
cluster = templateDict.get(cluster_key, 'auto').lower()
if cluster == 'auto':
cluster = templateAutoDict[cluster_key]
if cluster != 'local':
templateAutoDict['mintpy.compute.numWorker'] = '40'
## Update auto value of input template dict
for key, value in templateDict.items():
if value == 'auto' and key in templateAutoDict.keys():
templateDict[key] = templateAutoDict[key]
# Change yes --> True, no --> False and none --> None
special_values = {
'yes' : True,
'true' : True,
'no' : False,
'false': False,
'none' : None,
}
for key, value in templateDict.items():
value = value.lower()
if value in special_values.keys():
templateDict[key] = special_values[value]
return templateDict
def run_deramp(fname, ramp_type, mask_file=None, out_file=None, datasetName=None,
save_ramp_coeff=False, extra_meta=None):
""" Remove ramp from each 2D matrix of input file
Parameters: fname - str, data file to be deramped
ramp_type - str, name of ramp to be estimated.
mask_file - str, file of mask of pixels used for ramp estimation
out_file - str, output file name
datasetName - str, output dataset name, for ifgramStack file type only
save_ramp_coeff - bool, save the estimated ramp coefficients to text file
extra_meta - dict, extra metadata to add to the output file
Returns: out_file - str, output file name
"""
start_time = time.time()
# file/dir
fdir = os.path.dirname(fname)
fbase, fext = os.path.splitext(os.path.basename(fname))
# metadata
atr = readfile.read_attribute(fname)
ftype = atr['FILE_TYPE']
length = int(atr['LENGTH'])
width = int(atr['WIDTH'])
print(f'remove {ramp_type} ramp from file: {fname}')
out_file = out_file if out_file else os.path.join(fdir, f'{fbase}_ramp{fext}')
# ignore out_file for ifgramStack (write back to the same HDF5 file)
if ftype == 'ifgramStack':
out_file = fname
# mask
if os.path.isfile(mask_file):
mask = readfile.read(mask_file)[0]
print('read mask file: '+mask_file)
else:
mask = np.ones((length, width), dtype=np.bool_)
print('use mask of the whole area')
# write coefficient of specified surface function fit
coeff_file = None
if save_ramp_coeff:
coeff_file = os.path.join(fdir, f'rampCoeff_{fbase}.txt')
with open(coeff_file, 'w') as f:
f.write(f'# input file: {fname}\n')
f.write(f'# output file: {out_file}\n')
f.write(f'# ramp type: {ramp_type}\n')
# deramping
if ftype == 'timeseries':
# write HDF5 file with defined metadata and (empty) dataset structure
writefile.layout_hdf5(out_file, ref_file=fname, print_msg=True)
print('estimating phase ramp one date at a time ...')
date_list = timeseries(fname).get_date_list()
num_date = len(date_list)
prog_bar = ptime.progressBar(maxValue=num_date)
for i in range(num_date):
if coeff_file:
# prepend epoch name to line of coefficients
with open(coeff_file, 'a') as f:
f.write(f'{(date_list[i])} ')
# read
data = readfile.read(fname, datasetName=date_list[i])[0]
# deramp
data = deramp(
data,
mask,
ramp_type=ramp_type,
metadata=atr,
coeff_file=coeff_file,
)[0]
# write
writefile.write_hdf5_block(
out_file,
data=data,
datasetName='timeseries',
block=[i, i+1, 0, length, 0, width],
print_msg=False,
)
prog_bar.update(i+1, suffix=f'{i+1}/{num_date}')
prog_bar.close()
print(f'finished writing to file: {out_file}')
elif ftype == 'ifgramStack':
obj = ifgramStack(fname)
obj.open(print_msg=False)
if not datasetName:
datasetName = 'unwrapPhase'
with h5py.File(fname, 'a') as f:
ds = f[datasetName]
dsNameOut = f'{datasetName}_ramp'
if dsNameOut in f.keys():
dsOut = f[dsNameOut]
print(f'access HDF5 dataset /{dsNameOut}')
else:
dsOut = f.create_dataset(
dsNameOut,
shape=(obj.numIfgram, length, width),
dtype=np.float32,
chunks=True,
compression=None)
print(f'create HDF5 dataset /{dsNameOut}')
prog_bar = ptime.progressBar(maxValue=obj.numIfgram)
for i in range(obj.numIfgram):
if coeff_file:
# prepend IFG date12 to line of coefficients
with open(coeff_file, 'a') as f:
f.write(f'{str(obj.date12List[i])} ')
# read
data = ds[i, :, :]
# deramp
data = deramp(
data,
mask,
ramp_type=ramp_type,
metadata=atr,
coeff_file=coeff_file,
)[0]
# write
dsOut[i, :, :] = data
prog_bar.update(i+1, suffix=f'{i+1}/{obj.numIfgram}')
prog_bar.close()
print(f'finished writing to file: {fname}')
# Single Dataset File
else:
if coeff_file:
# prepend file-type to line of coefficients
with open(coeff_file, 'a') as f:
f.write('{} '.format(atr['FILE_TYPE']))
# read
if not datasetName and ftype == 'velocity':
datasetName = 'velocity'
data = readfile.read(fname, datasetName=datasetName)[0]
# deramp
data = deramp(
data,
mask,
ramp_type=ramp_type,
metadata=atr,
coeff_file=coeff_file,
)[0]
# write
print(f'writing >>> {out_file}')
writefile.write(data, out_file=out_file, ref_file=fname)
# add extra_meta to the output file
if extra_meta:
print('add/update the following metadata to file:')
add_attribute(out_file, extra_meta, print_msg=True)
# used time
m, s = divmod(time.time()-start_time, 60)
print(f'time used: {m:02.0f} mins {s:02.1f} secs.')
return out_file