Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ionex: migrate from interp2d to RegularGridInterpolator #1223

Merged
merged 1 commit into from
Jun 25, 2024
Merged

ionex: migrate from interp2d to RegularGridInterpolator #1223

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 20 additions & 22 deletions src/mintpy/objects/ionex.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
import re

import numpy as np
from scipy import interpolate
from scipy.interpolate import RegularGridInterpolator, interpn

from mintpy.utils import ptime

Expand Down Expand Up @@ -139,8 +139,8 @@ def interp_3d_rotate(interpfs, mins, lats, lons, utc_min, lat, lon):
ind1 = ind0 + 1
lon0 = lon + (utc_min - mins[ind0]) * 360. / (24. * 60.)
lon1 = lon + (utc_min - mins[ind1]) * 360. / (24. * 60.)
tec_val0 = interpfs[ind0](lon0, lat)
tec_val1 = interpfs[ind1](lon1, lat)
tec_val0 = interpfs[ind0]((lon0, lat))
tec_val1 = interpfs[ind1]((lon1, lat))
tec_val = ( (mins[ind1] - utc_min) / (mins[ind1] - mins[ind0]) * tec_val0
+ (utc_min - mins[ind0]) / (mins[ind1] - mins[ind0]) * tec_val1 )
return tec_val
Expand All @@ -160,28 +160,26 @@ def interp_3d_rotate(interpfs, mins, lats, lons, utc_min, lat, lon):
tec_val = np.zeros(num_pts, dtype=np.float32)
prog_bar = ptime.progressBar(maxValue=num_pts, print_msg=print_msg)
for i in range(num_pts):
tec_val[i] = interpolate.interp2d(
x=lons,
y=lats,
z=tec_maps[time_ind[i], :, :],
kind='linear',
)(lon[i], lat[i])
tec_val[i] = RegularGridInterpolator(
points=(lons, lats),
values=tec_maps[time_ind[i], :, :].T,
method='linear',
)((lon[i], lat[i]))

prog_bar.update(i+1, every=200)
prog_bar.close()

else:
tec_val = interpolate.interp2d(
x=lons,
y=lats,
z=tec_maps[time_ind[0], :, :],
kind='linear',
)(lon, lat)
tec_val = RegularGridInterpolator(
points=(lons, lats),
values=tec_maps[time_ind[0], :, :].T,
method='linear',
)((lon, lat))

elif interp_method in ['linear3d', 'trilinear']:
if not rotate_tec_map:
# option 1: interpolate between consecutive TEC maps
tec_val = interpolate.interpn(
tec_val = interpn(
points=(mins, np.flip(lats), lons),
values=np.flip(tec_maps, axis=1),
xi=(utc_min, lat, lon),
Expand All @@ -196,14 +194,14 @@ def interp_3d_rotate(interpfs, mins, lats, lons, utc_min, lat, lon):
interpfs = []
for i in range(len(mins)):
interpfs.append(
interpolate.interp2d(
x=lons,
y=lats,
z=tec_maps[i, :, :],
kind='linear',
RegularGridInterpolator(
points=(lons, lats),
values=tec_maps[i, :, :].T,
method='linear',
),
)

# interpolate
if isinstance(utc_min, np.ndarray):
num_pts = len(utc_min)
tec_val = np.zeros(num_pts, dtype=np.float32)
Expand All @@ -222,7 +220,7 @@ def interp_3d_rotate(interpfs, mins, lats, lons, utc_min, lat, lon):
interpfs,
mins, lats, lons,
utc_min, lat, lon,
)[0]
)

else:
msg = f'Un-recognized interp_method input: {interp_method}!'
Expand Down
4 changes: 4 additions & 0 deletions tests/objects/ionex.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ def prep_test_data(prep_mode=True):
def test_read_ionex():
"""Test mintpy.objects.ionex.read_ionex()"""

print('Test 1: read GIM files via ionex.read_ionex()')
time_ind = 1
x0, x1, y0, y1 = 3, 9, 28, 33
tec_aoi = np.array(
Expand All @@ -69,11 +70,13 @@ def test_read_ionex():

# compare
assert np.allclose(tec_maps[time_ind, y0:y1, x0:x1], tec_aoi)
print('Pass.')


def test_get_ionex_value():
"""Test mintpy.objects.ionex.get_ionex_value()"""

print('Test 2: Get the ionex for a given location/time via ionex.get_ionex_value()')
lat, lon = -21.3, -67.4 # northern Chile
utc_sec = 23 * 3600 + 7 * 60 # 23:07 UTC

Expand All @@ -96,6 +99,7 @@ def test_get_ionex_value():
rotate_tec_map=rotate,
)
assert np.allclose(tec_val, value, atol=1e-05, rtol=1e-05)
print(f'Interpolation method ({method:8s}) with rotation ({str(rotate):5s}): Pass.')


if __name__ == '__main__':
Expand Down