Closed Sugeradi985 closed 1 year ago
yunjunz
commented
1 year ago I would suggest check the following:
mintpy.troposphericDelay.method = gacos to turn on gacos tropospheric correction in mintpy.GACOS.h5, not ERA5.h5 file.
Sugeradi985
commented
1 year ago I would suggest check the following:
- you should set
mintpy.troposphericDelay.method = gacosto turn on gacos tropospheric correction in mintpy.- the gacos result is saved in
GACOS.h5, notERA5.h5file.
Hi, your response is greatly appreciated. So, you mean that the removal of atmospheric delay can only be done now using gacos. Am I correct? Considering how well I have used era-5 in the past.
Anyway, I'll give it a try using GACOS.
Thanks again.
yunjunz
commented
1 year ago So, you mean that the removal of atmospheric delay can only be done now using gacos. Am I correct?
No, all methods should work as described. I was having the impression that you were trying to use gacos for tropospheric correction.
While testing PyAPS/ERA5 on the Ridgecrest EQ example dataset, I had an error as well, which is now fixed by #1061. But the error I had is different from yours, could you test that PR to check if it fixes your issue as well?
Sugeradi985
commented
1 year ago So, you mean that the removal of atmospheric delay can only be done now using gacos. Am I correct?
No, all methods should work as described. I was having the impression that you were trying to use
gacosfor tropospheric correction.While testing PyAPS/ERA5 on the Ridgecrest EQ example dataset, I had an error as well, which is now fixed by #1061. But the error I had is different from yours, could you test that PR to check if it fixes your issue as well?
Hi, thanks for your prompt response. OK, I will try the Ridgecrest dataset then.
Sugeradi985
commented
1 year ago So, you mean that the removal of atmospheric delay can only be done now using gacos. Am I correct?
No, all methods should work as described. I was having the impression that you were trying to use
gacosfor tropospheric correction. While testing PyAPS/ERA5 on the Ridgecrest EQ example dataset, I had an error as well, which is now fixed by #1061. But the error I had is different from yours, could you test that PR to check if it fixes your issue as well?Hi, thanks for your prompt response. OK, I will try the Ridgecrest dataset then.
I attempted to correct the troposphere on the Ridgecrest dataset. I downloaded the grb files online and attempted to generate the ERA5.h5 file. But the same error message was received.
******************** step - correct_troposphere ********************
Input data seems to be geocoded. Lookup file not needed.
Atmospheric correction using Weather Re-analysis dataset (PyAPS, Jolivet et al., 2011)
Weather Re-analysis dataset: ERA5
tropo_pyaps3.py -f /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries.h5 --model ERA5 -g /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/inputs/geometryGeo.h5 -w /home/gmtsar-user/Share/process/WeatherData
weather model: ERA5 - dry (hydrostatic) and wet delay
weather directory: /home/gmtsar-user/Share/process/WeatherData
output tropospheric delay time-series file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/inputs/ERA5.h5
output corrected displacement time-series file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries_ERA5.h5
read dates/time info from file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries.h5
time of cloest available product: 14:00 UTC
--------------------------------------------------------------------------------
Download global atmospheric model files...
update mode: ON
output file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/inputs/ERA5.h5
1) output file either do NOT exist or is NOT newer than all GRIB files.
run or skip: run
--------------------------------------------------
downloading weather model data using PyAPS ...
common file size: 759240 bytes
number of grib files existed : 7
number of grib files to download: 0
--------------------------------------------------
--------------------------------------------------------------------------------
Calculate tropospheric delay and write to HDF5 file...
update mode: ON
output file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/inputs/ERA5.h5
1) output file either do NOT exist or is NOT newer than all GRIB files.
run or skip: run
open geometry file: geometryGeo.h5
reading incidenceAngle data from file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/inputs/geometryGeo.h5 ...
reading height data from file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/inputs/geometryGeo.h5 ...
--------------------------------------------------
create HDF5 file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/inputs/ERA5.h5 with w mode
create dataset : date of |S8 in size of (7,) with compression = None
create dataset : timeseries of <class 'numpy.float32'> in size of (7, 1125, 1500) with compression = None
close HDF5 file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/inputs/ERA5.h5
--------------------------------------------------
calculating absolute delay for each date using PyAPS (Jolivet et al., 2011; 2014) ...
number of grib files used: 7
Traceback (most recent call last):
File "/home/gmtsar-user/tools/miniconda3/envs/insar/bin/smallbaselineApp.py", line 8, in <module>
sys.exit(main())
File "/home/gmtsar-user/tools/MintPy/src/mintpy/cli/smallbaselineApp.py", line 208, in main
run_smallbaselineApp(inps)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/smallbaselineApp.py", line 1117, in run_smallbaselineApp
app.run(steps=inps.runSteps)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/smallbaselineApp.py", line 898, in run
self.run_tropospheric_delay_correction(sname)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/smallbaselineApp.py", line 648, in run_tropospheric_delay_correction
mintpy.cli.tropo_pyaps3.main(iargs)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/cli/tropo_pyaps3.py", line 166, in main
run_tropo_pyaps3(inps)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/tropo_pyaps3.py", line 675, in run_tropo_pyaps3
calc_delay_timeseries(inps)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/tropo_pyaps3.py", line 610, in calc_delay_timeseries
tropo_data = get_delay(
File "/home/gmtsar-user/tools/MintPy/src/mintpy/tropo_pyaps3.py", line 475, in get_delay
aps_obj = pa.PyAPS(
File "/home/gmtsar-user/tools/miniconda3/envs/insar/lib/python3.8/site-packages/pyaps3/objects.py", line 166, in __init__
hgt = np.linspace(self.dict['minAltP'], gph.max().round(), self.dict['nhgt'])
File "/home/gmtsar-user/tools/miniconda3/envs/insar/lib/python3.8/site-packages/numpy/core/_methods.py", line 40, in _amax
return umr_maximum(a, axis, None, out, keepdims, initial, where)
ValueError: zero-size array to reduction operation maximum which has no identity
Due to the fact that ERA5.h5 is only 10+kb, I believe there is an issue.
Then, I tested the ERA5.h5 file provided with the Ridgecrest dataset and it worked well.
******************** step - correct_troposphere ********************
Input data seems to be geocoded. Lookup file not needed.
Atmospheric correction using Weather Re-analysis dataset (PyAPS, Jolivet et al., 2011)
Weather Re-analysis dataset: ERA5
--------------------------------------------
Use existed tropospheric delay file: ./inputs/ERA5.h5
diff.py /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries.h5 ./inputs/ERA5.h5 -o /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries_ERA5.h5 --force
/home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries.h5 - ['./inputs/ERA5.h5'] --> /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries_ERA5.h5
the 1st input file is: timeseries
--------------------------------------------------
grab metadata from ref_file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries.h5
grab dataset structure from ref_file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries.h5
create HDF5 file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries_ERA5.h5 with w mode
create dataset : bperp of float32 in size of (7,) with compression = None
create dataset : date of |S8 in size of (7,) with compression = None
create dataset : timeseries of float32 in size of (7, 1125, 1500) with compression = None
close HDF5 file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries_ERA5.h5
read from file: ./inputs/ERA5.h5
* referencing data from ERA5.h5 to y/x: 1063/750
* referencing data from ERA5.h5 to date: 20190610
read from file: /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries.h5
--------------------------------------------------
open HDF5 file /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries_ERA5.h5 in a mode
writing dataset /timeseries block: [0, 7, 0, 1125, 0, 1500]
close HDF5 file /home/gmtsar-user/Documents/RidgecrestSenDT71/mintpy/timeseries_ERA5.h5.
time used: 00 mins 0.5 secsHmm... Could you check if your pyaps is installed and runs properly, by running the PyAPS/tests/test_calc.py script?
The output should look like below:
(insar) yunjunz:~>$ python tools/PyAPS/tests/test_calc.py
------------------------------------------------
import pyaps3 from /Users/yunjunz/tools/PyAPS/src/pyaps3/__init__.py
------------------------------------------------
test tropospheric delay calculation from ERA5.
------------------------------------------------
read ISCE geometry files: hgt/los/lat/lon.rdr
calculate tropospheric delay from GRB files...
------------------------------------------------
INFO: INCIDENCE ANGLE AS AN ARRAY
INFO: AREA COVERAGE IN SNWE: (33.85, 30.05, 129.05, 132.45)
PROGRESS: READING GRIB FILE
INFO: USING PRESSURE LEVELS OF ERA-INT OR ERA-5 DATA
INFO: IMAGE DIMENSIONS: 15 LATITUDES AND 13 LONGITUDES
PROGRESS: INTERPOLATING FROM PRESSURE TO HEIGHT LEVELS
PROGRESS: COMPUTING DELAY FUNCTIONS
INFO: INCIDENCE ANGLE AS AN ARRAY
INFO: AREA COVERAGE IN SNWE: (33.85, 30.05, 129.05, 132.45)
PROGRESS: READING GRIB FILE
INFO: USING PRESSURE LEVELS OF ERA-INT OR ERA-5 DATA
INFO: IMAGE DIMENSIONS: 15 LATITUDES AND 13 LONGITUDES
PROGRESS: INTERPOLATING FROM PRESSURE TO HEIGHT LEVELS
PROGRESS: COMPUTING DELAY FUNCTIONS
PROGRESS: FINE INTERPOLATION OF HEIGHT LEVELS
PROGRESS: CREATE THE BILINEAR INTERPOLATION FUNCTION
PROGRESS: MAPPING THE DELAY
[============================================================] 0s / 0s
PROGRESS: FINE INTERPOLATION OF HEIGHT LEVELS
PROGRESS: CREATE THE BILINEAR INTERPOLATION FUNCTION
PROGRESS: MAPPING THE DELAY
[============================================================] 0s / 0s
------------------------------------------------
Passed tropospheric delay calculation from ERA5.
------------------------------------------------Hmm... Could you check if your pyaps is installed and runs properly, by running the
PyAPS/tests/test_calc.pyscript?The output should look like below:
(insar) yunjunz:~>$ python tools/PyAPS/tests/test_calc.py ------------------------------------------------ import pyaps3 from /Users/yunjunz/tools/PyAPS/src/pyaps3/__init__.py ------------------------------------------------ test tropospheric delay calculation from ERA5. ------------------------------------------------ read ISCE geometry files: hgt/los/lat/lon.rdr calculate tropospheric delay from GRB files... ------------------------------------------------ INFO: INCIDENCE ANGLE AS AN ARRAY INFO: AREA COVERAGE IN SNWE: (33.85, 30.05, 129.05, 132.45) PROGRESS: READING GRIB FILE INFO: USING PRESSURE LEVELS OF ERA-INT OR ERA-5 DATA INFO: IMAGE DIMENSIONS: 15 LATITUDES AND 13 LONGITUDES PROGRESS: INTERPOLATING FROM PRESSURE TO HEIGHT LEVELS PROGRESS: COMPUTING DELAY FUNCTIONS INFO: INCIDENCE ANGLE AS AN ARRAY INFO: AREA COVERAGE IN SNWE: (33.85, 30.05, 129.05, 132.45) PROGRESS: READING GRIB FILE INFO: USING PRESSURE LEVELS OF ERA-INT OR ERA-5 DATA INFO: IMAGE DIMENSIONS: 15 LATITUDES AND 13 LONGITUDES PROGRESS: INTERPOLATING FROM PRESSURE TO HEIGHT LEVELS PROGRESS: COMPUTING DELAY FUNCTIONS PROGRESS: FINE INTERPOLATION OF HEIGHT LEVELS PROGRESS: CREATE THE BILINEAR INTERPOLATION FUNCTION PROGRESS: MAPPING THE DELAY [============================================================] 0s / 0s PROGRESS: FINE INTERPOLATION OF HEIGHT LEVELS PROGRESS: CREATE THE BILINEAR INTERPOLATION FUNCTION PROGRESS: MAPPING THE DELAY [============================================================] 0s / 0s ------------------------------------------------ Passed tropospheric delay calculation from ERA5. ------------------------------------------------
After thoroughly examining the detailed error message, I believe it was not caused by mintpy or pyaps, but by an environmental configuration issue.
I believe that numpy or scipy, or something else, is the cause of the issue, although it's not clear.
So I re-configured the Conda environment and it now works well.
Sorry for my own reasons, I mistakenly thought it was a bug.
I believe we can close this comment this time.
Thank you very much!
yunjunz
commented
1 year ago No worries, I am glad to hear you locate the cause and fix the issue!
Sugeradi985
commented
1 year ago No worries, I am glad to hear you locate the cause and fix the issue!
Hi, it seems that the error is still present.
Traceback (most recent call last):
File "/home/gmtsar-user/tools/mambaforge/envs/insar/bin/smallbaselineApp.py", line 8, in <module>
sys.exit(main())
File "/home/gmtsar-user/tools/MintPy/src/mintpy/cli/smallbaselineApp.py", line 208, in main
run_smallbaselineApp(inps)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/smallbaselineApp.py", line 1117, in run_smallbaselineApp
app.run(steps=inps.runSteps)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/smallbaselineApp.py", line 898, in run
self.run_tropospheric_delay_correction(sname)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/smallbaselineApp.py", line 648, in run_tropospheric_delay_correction
mintpy.cli.tropo_pyaps3.main(iargs)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/cli/tropo_pyaps3.py", line 166, in main
run_tropo_pyaps3(inps)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/tropo_pyaps3.py", line 675, in run_tropo_pyaps3
calc_delay_timeseries(inps)
File "/home/gmtsar-user/tools/MintPy/src/mintpy/tropo_pyaps3.py", line 610, in calc_delay_timeseries
tropo_data = get_delay(
File "/home/gmtsar-user/tools/MintPy/src/mintpy/tropo_pyaps3.py", line 475, in get_delay
aps_obj = pa.PyAPS(
File "/home/gmtsar-user/tools/PyAPS/src/pyaps3/objects.py", line 177, in __init__
hgt = np.linspace(self.dict['minAltP'], gph.max().round(), self.dict['nhgt'])
File "/home/gmtsar-user/tools/mambaforge/envs/insar/lib/python3.9/site-packages/numpy/core/_methods.py", line 41, in _amax
return umr_maximum(a, axis, None, out, keepdims, initial, where)
ValueError: zero-size array to reduction operation maximum which has no identityI'm trying to print the message from some scripts.
PyAPS/src/pyaps3/objects.py", line 177, in init hgt = np.linspace(self.dict['minAltP'], gph.max().round(), self.dict['nhgt'])
Referring to the error message "ValueError: zero-size array to reduction operation maximum which has no identity", I think it was caused by the 'gph.max().round' function.
Then I realized that the gph (array) was empty, and defined in era.py line 208:
lats, lons = grb.latlons()
#if model == 'ERA5':
# lons[lons < 0.] += 360.
g = cdic['g']
#extract indices
mask = ((lats > minlat) & (lats < maxlat)) & ((lons > minlon) & (lons < maxlon))
uu = [i for i in list(range(np.shape(mask)[0])) if any(mask[i,:])]
vv = [j for j in list(range(np.shape(mask)[1])) if any(mask[:,j])]
latlist = lats[uu,:][:,vv]
lonlist = lons[uu,:][:,vv]
nlat, nlon = latlist.shape
####Create arrays for 3D storage
gph = np.zeros((nlvls, nlat, nlon)) #Potential heightAfter following the output message step by step, I discovered that the coverage of lats(array) is different with the true latitude of the dataset.
With the Ridgecrest dataset, the output message is as follows:
calculating absolute delay for each date using PyAPS (Jolivet et al., 2011; 2014) ...
number of grib files used: 7
lats[10]=
[37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5
37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5
37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5
37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5
37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5
37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5]
lons[10]=
[-130. -129.75 -129.5 -129.25 -129. -128.75 -128.5 -128.25 -128.
-127.75 -127.5 -127.25 -127. -126.75 -126.5 -126.25 -126. -125.75
-125.5 -125.25 -125. -124.75 -124.5 -124.25 -124. -123.75 -123.5
-123.25 -123. -122.75 -122.5 -122.25 -122. -121.75 -121.5 -121.25
-121. -120.75 -120.5 -120.25 -120. -119.75 -119.5 -119.25 -119.
-118.75 -118.5 -118.25 -118. -117.75 -117.5 -117.25 -117. -116.75
-116.5 -116.25 -116. -115.75 -115.5 -115.25 -115. -114.75 -114.5
-114.25 -114. -113.75 -113.5 -113.25 -113. -112.75 -112.5 -112.25
-112. -111.75 -111.5 -111.25 -111. -110.75 -110.5 -110.25 -110. ]
mask:
[[False False False ... False False False]
[False False False ... False False False]
[False False False ... False False False]
...
[False False False ... False False False]
[False False False ... False False False]
[False False False ... False False False]]
minlat,minlon,maxlat,maxlon:
-1.199681364756279 , -122.68980303959215 , 1.2003260023253484 , -120.28979109842726
nlvls:
37
nlat:
0
nlon:
0
uu and vv
[] , []We can see that, the output minlat (-1.199681364756279), maxlat(1.2003260023253484 ) are incorrect. So the mask (array) is full of False.
These two arrays are calculated using the formula in PyAPS/src/pyaps3/objects.py, in lin 111:
# Problems in isce when lon and lat arrays have weird numbers
if self.grib in ('ERA5','ERAINT','HRES'):
self.lon[self.lon > 180.] -= 360.
else:
self.lon[self.lon < 0.] += 360.0
self.minlon = np.nanmin(self.lon[np.nonzero(self.mask)]) - self.bufspc
self.maxlon = np.nanmax(self.lon[np.nonzero(self.mask)]) + self.bufspc
self.minlat = np.nanmin(self.lat[np.nonzero(self.mask)]) - self.bufspc
self.maxlat = np.nanmax(self.lat[np.nonzero(self.mask)]) + self.bufspc
if verb:
print('INFO: AREA COVERAGE IN SNWE: ({:.2f}, {:.2f}, {:.2f}, {:.2f})'.format(
self.maxlat, self.minlat, self.minlon, self.maxlon))So, I think there may be some errors in calculating the min- and max-latitude.
yunjunz
commented
1 year ago For errors in pyaps, please consider opening an issue on its repo instead. Nevertheless, I could successfully run Ridgecrest dataset with pyaps, so my guess is that your code is not properly installed. I would recommend re-installing them from scratch before digging into the code. Both PyAPS and MintPy have tests, use them to check your installation.
Sugeradi985
commented
1 year ago Hmm, I reinstalled them yesterday following your guide on conda-envs. Both tests are effective.
May I ask you to test the Ridgecrest dataset without the ERA5.h5 file under mintpy/inputs ? Simply delete the mintpy/inputs/ERA5.h5 fisrt, then run pyaps.
Thanks!
yunjunz
commented
1 year ago While testing PyAPS/ERA5 on the Ridgecrest EQ example dataset, I had an error as well, which is now fixed by https://github.com/insarlab/MintPy/pull/1061.
I already did previously. If you could strictly follow the bug report template, and fill in all the details, it could help to replicate your issue.
Problem description
Hi, wish you all the best.
I want to do the time series analysis on the iterferograms generated by Hyp3.
And I have pre-processed the tiffs accordingly.
When came to the "correct_troposphere" ,I encounter errors.
After checking the generated file "ERA5.h5," which was just 14.4 kilobytes, it became apparent that something was wrong.
I used info.py to view the information for "ERA5.h5", but it was full of zeros.
I used view.py to plot the "ERA5.h5", but it was empty.
I've tried using different sets of interferograms (also generated by Hyp3), and they all come out the same way.
Full script that generated the error
Full error message
System information
Any help is greatly appreciated.