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.

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

Change the way observer location is calculated in Map #3115

Merged
merged 4 commits into from
May 21, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
4 changes: 4 additions & 0 deletions changelog/3115.breaking.1.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
`sunpy.map.GenericMap` now checks for a complete observer location rather than
individually defaulting coordinates (lat, lon, distance) to Earth position. If
any one of the three coordinates is missing from the header the observer will
be defaulted to Earth and a warning raised.
1 change: 1 addition & 0 deletions changelog/3115.breaking.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
`sunpy.map.GenericMap` will no longer use the key `solar_b0` as a value for heliographic latitude.
121 changes: 56 additions & 65 deletions sunpy/map/mapbase.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@
from sunpy.image.transform import affine_transform
from sunpy.image.resample import reshape_image_to_4d_superpixel
from sunpy.image.resample import resample as sunpy_image_resample
from sunpy.coordinates import get_sun_B0, get_sun_L0, get_sunearth_distance
from sunpy.coordinates import get_earth
from sunpy.util import expand_list
from sunpy.util.exceptions import SunpyUserWarning

from astropy.nddata import NDData
Expand Down Expand Up @@ -421,20 +422,6 @@ def detector(self):
"""Detector name."""
return self.meta.get('detector', "")

@property
def dsun(self):
"""The observer distance from the Sun."""
dsun = self.meta.get('dsun_obs', None)

if dsun is None:
if self._default_dsun is None:
warnings.warn("Missing metadata for Sun-spacecraft separation: assuming Sun-Earth distance.",
SunpyUserWarning)
self._default_dsun = get_sunearth_distance(self.date).to(u.m)
return self._default_dsun

return u.Quantity(dsun, 'm')

@property
def exposure_time(self):
"""Exposure time of the image in seconds."""
Expand Down Expand Up @@ -562,71 +549,75 @@ def coordinate_system(self):
self.meta.get('ctype2', 'HPLT- '))

@property
def carrington_longitude(self):
"""Carrington longitude (crln_obs)."""
carrington_longitude = self.meta.get('crln_obs', None)

if carrington_longitude is None:
if self._default_carrington_longitude is None:
warnings.warn("Missing metadata for Carrington longitude: assuming Earth-based observer.",
SunpyUserWarning)
self._default_carrington_longitude = get_sun_L0(self.date)
carrington_longitude = self._default_carrington_longitude
def _supported_observer_coordinates(self):
"""
A list of supported coordinate systems.

if isinstance(carrington_longitude, str):
carrington_longitude = float(carrington_longitude)
This is a list so it can easily maintain a strict order. The list of
two element tuples, the first item in the tuple is the keys that need
to be in the header to use this coordinate system and the second is the
kwargs to SkyCoord.
"""
return [(('hgln_obs', 'hglt_obs', 'dsun_obs'), {'lon': self.meta.get('hgln_obs'),
'lat': self.meta.get('hglt_obs'),
'radius': self.meta.get('dsun_obs'),
'unit': (u.deg, u.deg, u.m),
'frame': "heliographic_stonyhurst"}),
(('crln_obs', 'crlt_obs', 'dsun_obs'), {'lon': self.meta.get('crln_obs'),
'lat': self.meta.get('crlt_obs'),
'radius': self.meta.get('dsun_obs'),
'unit': (u.deg, u.deg, u.m),
'frame': "heliographic_carrington"}),]

return u.Quantity(carrington_longitude, 'deg')
def _remove_existing_observer_location(self):
"""
Remove all keys that this map might use for observer location.
"""
all_keys = expand_list([e[0] for e in self._supported_observer_coordinates])
for key in all_keys:
self.meta.pop(key)

@property
def heliographic_latitude(self):
"""Heliographic latitude."""
heliographic_latitude = self.meta.get('hglt_obs',
self.meta.get('crlt_obs',
self.meta.get('solar_b0', None)))
def observer_coordinate(self):
"""
The Heliographic Stonyhurst Coordinate of the observer.
"""
for keys, kwargs in self._supported_observer_coordinates:
if all([k in self.meta for k in keys]):
return SkyCoord(obstime=self.date, **kwargs).heliographic_stonyhurst

if heliographic_latitude is None:
if self._default_heliographic_latitude is None:
warnings.warn("Missing metadata for heliographic latitude: assuming Earth-based observer.",
SunpyUserWarning)
self._default_heliographic_latitude = get_sun_B0(self.date)
heliographic_latitude = self._default_heliographic_latitude
all_keys = [str(e[0]) for e in self._supported_observer_coordinates]
all_keys = '\n'.join(all_keys)
warning_message = ("Missing metadata for observer: assuming Earth-based observer."
Copy link
Member

Choose a reason for hiding this comment

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

I think it would be useful to tell the user exactly which metadata is missing here (instead of just telling them that something is missing).

Copy link
Member Author

Choose a reason for hiding this comment

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

It's not trivial to work out what ones are missing, i.e. you could have a weird header with a mix of crln and hglt and it would be hard to work out which one the user was aiming for.

I will give it some thought.

Copy link
Member

Choose a reason for hiding this comment

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

I think improved docs on map meta will really help here

Copy link
Member Author

Choose a reason for hiding this comment

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

I opened an issue to keep track of this because I didn't want to block this on it.

"The following sets of keys were checked:\n" + all_keys)
warnings.warn(warning_message, SunpyUserWarning)

if isinstance(heliographic_latitude, str):
heliographic_latitude = float(heliographic_latitude)
return get_earth(self.date)

return u.Quantity(heliographic_latitude, 'deg')
@property
def heliographic_latitude(self):
"""Heliographic latitude."""
return self.observer_coordinate.lat

@property
def heliographic_longitude(self):
"""Heliographic longitude."""
heliographic_longitude = self.meta.get('hgln_obs', None)

if heliographic_longitude is None:
if self.meta.get('crln_obs', None) is not None:
heliographic_longitude = self.meta['crln_obs'] * u.deg - get_sun_L0(self.date)
else:
if self._default_heliographic_longitude is None:
warnings.warn("Missing metadata for heliographic longitude: assuming longitude of 0 degrees.",
SunpyUserWarning)
self._default_heliographic_longitude = 0
heliographic_longitude = self._default_heliographic_longitude
return self.observer_coordinate.lon

if isinstance(heliographic_longitude, str):
heliographic_longitude = float(heliographic_longitude)
@property
def carrington_latitude(self):
"""Carrington latitude."""
return self.observer_coordinate.heliographic_carrington.lat

return u.Quantity(heliographic_longitude, 'deg')
@property
def carrington_longitude(self):
"""Carrington longitude."""
return self.observer_coordinate.heliographic_carrington.lon

@property
def observer_coordinate(self):
"""
The Heliographic Stonyhurst Coordinate of the observer.
"""
return SkyCoord(lat=self.heliographic_latitude,
lon=self.heliographic_longitude,
radius=self.dsun,
obstime=self.date,
frame='heliographic_stonyhurst')
def dsun(self):
"""The observer distance from the Sun."""
return self.observer_coordinate.radius.to('m')

@property
def _reference_longitude(self):
Expand Down
45 changes: 8 additions & 37 deletions sunpy/map/sources/sdo.py
Original file line number Diff line number Diff line change
Expand Up @@ -63,43 +63,14 @@ def __init__(self, data, header, **kwargs):
self.plot_settings['norm'] = ImageNormalize(stretch=source_stretch(self.meta, AsinhStretch(0.01)))

@property
def observer_coordinate(self):
"""
The Heliographic Stonyhurst Coordinate of the observer.

This coordinate is determined using the Heliocentric Aries Ecliptic (HAE) coordinates
in the header.
"""
vector = CartesianRepresentation(self.meta['haex_obs'],
self.meta['haey_obs'],
self.meta['haez_obs'])
coord = SkyCoord(vector * u.m, frame=HeliocentricMeanEcliptic, obstime=self.date)
return coord.heliographic_stonyhurst

@property
def heliographic_latitude(self):
"""Heliographic latitude."""
return self.observer_coordinate.lat

@property
def heliographic_longitude(self):
"""Heliographic longitude."""
return self.observer_coordinate.lon

@property
def carrington_latitude(self):
"""Carrington latitude."""
return self.observer_coordinate.heliographic_carrington.lat

@property
def carrington_longitude(self):
"""Carrington longitude."""
return self.observer_coordinate.heliographic_carrington.lon

@property
def dsun(self):
"""The observer distance from the Sun."""
return self.observer_coordinate.radius.to('m')
def _supported_observer_coordinates(self):
return [(('haex_obs', 'haey_obs', 'haez_obs'), {'x': self.meta.get('haex_obs'),
'y': self.meta.get('haey_obs'),
'z': self.meta.get('haez_obs'),
'unit': u.m,
'representation': CartesianRepresentation,
'frame': HeliocentricMeanEcliptic})
] + super()._supported_observer_coordinates

@property
def observatory(self):
Expand Down
19 changes: 13 additions & 6 deletions sunpy/map/tests/test_mapbase.py
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ def test_detector(generic_map):


def test_dsun(generic_map):
with pytest.warns(SunpyUserWarning, match='Missing metadata for Sun-spacecraft separation'):
with pytest.warns(SunpyUserWarning, match='Missing metadata for observer: assuming Earth-based observer.*'):
assert generic_map.dsun == sunpy.coordinates.get_sunearth_distance(generic_map.date).to(u.m)


Expand All @@ -180,17 +180,17 @@ def test_coordinate_system(generic_map):


def test_carrington_longitude(generic_map):
with pytest.warns(SunpyUserWarning, match='Missing metadata for Carrington longitude'):
with pytest.warns(SunpyUserWarning, match='Missing metadata for observer: assuming Earth-based observer.*'):
assert generic_map.carrington_longitude == sunpy.coordinates.get_sun_L0(generic_map.date)


def test_heliographic_latitude(generic_map):
with pytest.warns(SunpyUserWarning, match='Missing metadata for heliographic latitude'):
with pytest.warns(SunpyUserWarning, match='Missing metadata for observer: assuming Earth-based observer.*'):
assert generic_map.heliographic_latitude == sunpy.coordinates.get_sun_B0(generic_map.date)


def test_heliographic_longitude(generic_map):
with pytest.warns(SunpyUserWarning, match='Missing metadata for heliographic longitude'):
with pytest.warns(SunpyUserWarning, match='Missing metadata for observer: assuming Earth-based observer.*'):
assert generic_map.heliographic_longitude == 0.


Expand All @@ -211,6 +211,13 @@ def test_heliographic_longitude_crln(hmi_test_map):
assert hmi_test_map.heliographic_longitude == hmi_test_map.carrington_longitude - sunpy.coordinates.get_sun_L0(hmi_test_map.date)


def test_remove_observers(aia171_test_map):
aia171_test_map._remove_existing_observer_location()
with pytest.warns(SunpyUserWarning,
match='Missing metadata for observer: assuming Earth-based observer.*'):
aia171_test_map.observer_coordinate


# ==============================================================================
# Test Rotation WCS conversion
# ==============================================================================
Expand Down Expand Up @@ -697,8 +704,8 @@ def test_missing_metadata_warnings():
header['cunit2'] = 'arcsec'
array_map = sunpy.map.Map(np.random.rand(20, 15), header)
array_map.peek()
# There should be 4 warnings for missing metadata
assert len([w for w in record if 'Missing metadata' in str(w)]) == 4
# There should be 2 warnings for missing metadata (obstime and observer location)
assert len([w for w in record if w.category is SunpyUserWarning]) == 2


def test_fits_header(aia171_test_map):
Expand Down