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

Fix get_observer_look for satellite exactly at nadir #77

Merged
merged 9 commits into from
Sep 23, 2021
Merged

Fix get_observer_look for satellite exactly at nadir #77

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
1efe676
Fix get_observer_look for satellite exactly at nadir.
Apr 12, 2021
1bf54ef
Handle problematic az_data for zero divided by zero
May 20, 2021
8023923
Extend test of satellite azimuth differense angle
May 21, 2021
89beb6d
Simplyfy get_observer_look
May 21, 2021
1efbd75
Change undefined azimuth value to 180
May 21, 2021
5e71510
Fix problem for macos (hopefully)
May 21, 2021
a885dce
Fix deepcode issue and reduce number of elements for test
May 21, 2021
7099f81
Attempt two to fix deepcode issue with random numbers
May 24, 2021
ca1d4d2
use the random number generator for all calls (forgot one)
May 24, 2021
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
30 changes: 26 additions & 4 deletions pyorbital/orbital.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,21 +133,43 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt):
else:
az_data = az_

# When elevation is 90 degrees both top_s and top_e are close to zero
# and azimuth angle is undefined.
undet_azi = np.logical_and(top_s == 0, top_e == 0)
if has_dask and isinstance(az_data, da.Array):
az_data = da.where(undet_azi, 0, az_data)
az_data = da.where(top_s > 0, az_data + np.pi, az_data)
az_data = da.where(az_data < 0, az_data + 2 * np.pi, az_data)
else:
az_data[np.where(undet_azi)] = 0
az_data[np.where(top_s > 0)] += np.pi
az_data[np.where(az_data < 0)] += 2 * np.pi
ninahakansson marked this conversation as resolved.
Show resolved Hide resolved

rg_ = np.sqrt(rx * rx + ry * ry + rz * rz)

top_z_divided_by_rg_ = top_z / rg_
el_ = np.arcsin(top_z_divided_by_rg_)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using clip might help avoid the extra filtering later:

Suggested change
top_z_divided_by_rg_ = top_z / rg_
el_ = np.arcsin(top_z_divided_by_rg_)
top_z_divided_by_rg_ = top_z / rg_
top_z_divided_by_rg_ = top_z_divided_by_rg_.clip(max=1)
el_ = np.arcsin(top_z_divided_by_rg_)


if has_xarray and isinstance(az_, xr.DataArray):
az_.data = az_data
el_data = el_.data
else:
az_ = az_data
el_data = el_

rg_ = np.sqrt(rx * rx + ry * ry + rz * rz)
el_ = np.arcsin(top_z / rg_)
# Due to rounding top_z can be larger than rg_ (when el_ ~ 90).
# And azimuth undefined when elevation is 90 degrees
if has_dask and isinstance(az_data, da.Array):
el_data = da.where(top_z_divided_by_rg_ > 1.0, np.pi/2, el_data)
az_data = da.where(el_data == np.pi / 2, 0, az_data)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of this filtering and the undet_azi higher up, would np.nan_to_num be helpfull here? or do we still want some azimuths to be nans?

else:
el_data[np.where(top_z_divided_by_rg_ > 1.0)] = np.pi/2
az_data[np.where(el_data == np.pi/2)] = 0
ninahakansson marked this conversation as resolved.
Show resolved Hide resolved

if has_xarray and isinstance(az_, xr.DataArray):
az_.data = az_data
el_.data = el_data
else:
el_ = el_data
az_ = az_data
return np.rad2deg(az_), np.rad2deg(el_)


Expand Down
8 changes: 4 additions & 4 deletions pyorbital/tests/test_orbital.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -220,15 +220,15 @@ class TestGetObserverLook(unittest.TestCase):
def setUp(self):
self.t = datetime(2018, 1, 1, 0, 0, 0)
self.sat_lon = np.array([[-89.5, -89.4], [-89.3, -89.2]])
self.sat_lat = np.array([[45.5, 45.4], [45.3, 45.2]])
self.sat_lat = np.array([[45.5, 45.4], [45.3, 40.2]])
self.sat_alt = np.array([[35786, 35786], [35786, 35786]])
self.lon = np.array([[-85.5, -85.4], [-85.3, -85.2]])
self.lon = np.array([[-85.5, -85.4], [-85.3, -89.2]])
self.lat = np.array([[40.5, 40.4], [40.3, 40.2]])
self.alt = np.zeros((2, 2))
self.exp_azi = np.array([[331.00275902, 330.95954165],
[330.91642994, 330.87342384]])
[330.91642994, 0]])
self.exp_elev = np.array([[83.18070976, 83.17788976],
[83.17507167, 83.1722555]])
[83.17507167, 90]])

def test_basic_numpy(self):
"""Test with numpy array inputs"""
Expand Down