Skip to content

map reprojection but better #38

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

Open
elastufka wants to merge 1 commit into
base: main
Choose a base branch
from
Open

map reprojection but better #38

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
136 changes: 94 additions & 42 deletions stixpy/vis/map_reprojection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,7 @@

from pathlib import Path

import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
import stixcore.data.test
Expand All @@ -87,8 +88,10 @@
from stixcore.data.test import test_data
from sunpy.coordinates import frames
from sunpy.coordinates import get_body_heliographic_stonyhurst
import drms
import warnings

__all__ = ['get_solo_position', 'reproject_map', 'create_headers', 'plot_map_reproj']
__all__ = ['get_solo_position', 'get_sdo_position', 'reproject_map', 'create_header', 'plot_map_reproj']


def get_solo_position(map):
Expand All @@ -103,63 +106,109 @@ def get_solo_position(map):
Returns
-------
solo_hgs : `astropy.coordinates.SkyCoord`
The position of SOLO in HeliographicStonyhurst frame
Reference coordinate in a helioprojective frame, with observer at position of SOLO
"""
p = Spice.instance.get_position(date=map.date.datetime, frame='SOLO_HEE')
solo_hee = SkyCoord(*p, frame=frames.HeliocentricEarthEcliptic, representation_type='cartesian',
obstime=map.date.datetime)
# Converting HeliocentricEarthEcliptic coords of SOLAR ORBITER position to
# HeliographicStonyhurst frame
solo_hgs = solo_hee.transform_to(frames.HeliographicStonyhurst)
return solo_hgs
return SkyCoord(0*u.arcsec,0*u.arcsec, obstime=map.date.datetime,observer=solo_hgs,frame='helioprojective')

def get_sdo_position(map):
"""
Return the position of SDO at the time the map was observed

Parameters
----------
map : `sunpy.map.Map`
Any solar map

Returns
-------
sdo_hgs : `astropy.coordinates.SkyCoord`
Reference coordinate in a helioprojective frame, with observer at position of SDO
"""
date_obs_str=map.meta['date-obs']
client = drms.Client()
kstr='CUNIT1,CUNIT2,CRVAL1,CRVAL2,CDELT1,CDELT2,CRPIX1,CRPIX2,CTYPE1,CTYPE2,HGLN_OBS,HGLT_OBS,DSUN_OBS,RSUN_OBS,DATE-OBS,RSUN_REF,CRLN_OBS,CRLT_OBS,EXPTIME,INSTRUME,WAVELNTH,WAVEUNIT,TELESCOP,LVL_NUM,CROTA2'
df = client.query(f'aia.lev1_uv_24s[{date_obs_str}/1m@30s]', key=kstr)
if df.empty: # Event not in database yet or other error
observer=get_observer(date_obs,obs='Earth',sc=sc)
warnings.warn(f"FITS headers for {date_obs} not available, using Earth observer" )
else:
try:
meta=df.where(df.WAVELNTH==1600).dropna().iloc[0].to_dict()
except IndexError:
meta=df.iloc[0].to_dict()
if np.isnan(meta['CRPIX1']):
meta['CRPIX1']=0.
if np.isnan(meta['CRPIX2']):
meta['CRPIX2']=0.
sdo_hgs = SkyCoord(meta['HGLN_OBS']*u.deg, meta['HGLT_OBS']*u.deg,meta['DSUN_OBS']*u.m, frame=frames.HeliographicStonyhurst, obstime=meta['DATE-OBS'])

return SkyCoord(0*u.arcsec,0*u.arcsec, obstime=map.date.datetime,observer=sdo_hgs,frame='helioprojective')


def create_headers(obs_ref_coord, map, out_shape=None, out_scale=None):
def create_header(obs_ref_coord, map, out_shape=None, out_scale=None):
"""
Generates MetaDict and WCS headers for reprojected map
Generates MetaDict header for reprojected map

Parameters
----------
obs_ref_coord : `sunpy.map.Map`
obs_ref_coord : `astropy.coordinates.SkyCoord`
Target WCS reference coordinate (as seen by observer).
Generated in map_reproject_to_observer function.
Generated in get_SPACECRAFT_position function.

map : `sunpy.map.Map`
Map to be reprojected.

out_shape : `tuple`
The shape of the reprojected map - defaults to (512, 512)
The shape of the reprojected map - defaults to input shape for full-disk input maps, otherwise will be calculated based on reprojected shape of submap

out_scale : `Quantity`
The scale of the output map
The scale of the output map. Defaults to 1":1px

Returns
-------
obs_wcs_header : `astropy.wcs`
WCS header for reprojected map
obs_metadict_header : `MetaDict`
MetaDict header for reprojected map
"""
if out_scale is None:
out_scale = u.Quantity(map.scale)

if out_shape is None:
out_shape = map.data.shape

if map.wavelength.unit.to_string() == '':
obs_metadict_header = sunpy.map.make_fitswcs_header(
out_shape,
obs_ref_coord,
scale=out_scale,
rotation_matrix=map.rotation_matrix,
instrument=map.detector)
else:
obs_metadict_header = sunpy.map.make_fitswcs_header(
out_shape,
obs_ref_coord,
scale=out_scale,
rotation_matrix=map.rotation_matrix,
instrument=map.detector,
wavelength=map.wavelength)
obs_wcs_header = WCS(obs_metadict_header)
return obs_wcs_header, obs_metadict_header
out_shape = (1,1) # So that shape can be modified later

target_observer = obs_ref_coord
target_wcs=WCS(sunpy.map.make_fitswcs_header(out_shape,target_observer)) #wcs_utils.reference_coordinate_from_frame(wcs_utils.solar_wcs_frame_mapping(target_wcs)) # Must be SkyCoord ie reference coordinate
edges_pix = np.concatenate(sunpy.map.map_edges(map))
edges_coord = map.pixel_to_world(edges_pix[:, 0], edges_pix[:, 1])
new_edges_coord = edges_coord.transform_to(target_observer)
new_edges_xpix, new_edges_ypix = target_wcs.world_to_pixel(new_edges_coord)

# Check for NaNs
if new_edges_xpix[np.isnan(new_edges_xpix)] != np.array([]) or new_edges_ypix[np.isnan(new_edges_ypix)] != np.array([]):
cc = sunpy.map.all_coordinates_from_map(map).transform_to(target_observer)
on_disk = sunpy.map.coordinate_is_on_solar_disk(cc)
on_disk_coordinates = cc[on_disk]
nan_percent = 1. - len(on_disk_coordinates)/len(cc.flatten())
if nan_percent > 0.5:
warnings.warn(f"{nan_percent*100:.1f}% of pixels in reprojected map are NaN!")

# Determine the extent needed - use of nanmax/nanmin means only on-disk coords are considered
left, right = np.nanmin(new_edges_xpix), np.nanmax(new_edges_xpix)
bottom, top = np.nanmin(new_edges_ypix), np.nanmax(new_edges_ypix)

# Adjust the CRPIX and NAXIS values
target_header = sunpy.map.make_fitswcs_header(out_shape, target_observer,instrument=map.detector) # What to do about scale?
target_header['crpix1'] -= left
target_header['crpix2'] -= bottom
target_header['naxis1'] = int(np.ceil(right - left))
target_header['naxis2'] = int(np.ceil(top - bottom))

return target_header


def reproject_map(map, observer, out_shape=None):
Expand All @@ -171,24 +220,27 @@ def reproject_map(map, observer, out_shape=None):
map : `sunpy.map.Map`
The input map to be reprojected
observer : `astropy.coordinates.SkyCoord`
The coordinates of the observer in HeliographicStonyhurst frame
The reference coordinate of the observer in HeliographicStonyhurst frame
out_shape : `tuple`
The shape of the reprojected map - defaults to same size as input map
The shape of the reprojected map - defaults to same size as input map if map is full-disk, otherwise the reprojected shape will be calculated based on the submap coordinates

Returns
-------
map : `sunpy.map.Map`
outmap : `sunpy.map.Map`
Reprojected map
"""
if out_shape is None:
if out_shape is None and sunpy.map.contains_full_disk(map):
out_shape = map.data.shape

obs_ref_coord = SkyCoord(map.reference_coordinate.Tx, map.reference_coordinate.Ty,
obstime=map.reference_coordinate.obstime,
observer=observer, frame="helioprojective")
obs_wcs_header, obs_metadict_header = create_headers(obs_ref_coord, map, out_shape=out_shape)
output, footprint = reproject_interp(map, obs_wcs_header, out_shape)
outmap = sunpy.map.Map(output, obs_metadict_header)

# Make sure observer is a reference coordinate, not just an observer coordinate
if not isinstance(observer, SkyCoord):
observer = SkyCoord(0*u.arcsec,0*u.arcsec, obstime=map.date.datetime,observer=observer,frame='helioprojective') # What if user inputs WCS? I know I wrote a observer_from_wcs somewhere ... find it and add here

obs_metadict_header = create_header(observer, map, out_shape=out_shape)
# Don't try to reproject to a huge array
if obs_metadict_header['naxis1'] > 2048 or obs_metadict_header['naxis2'] > 2048: # Reprojection will take at least 30s in this case.
warnings.warn(f"Reprojection to {obs_metadict_header['naxis1']}x{obs_metadict_header['naxis2']} array not permitted! Check that your desired observer is correct, use a submap of the region, or scale down. Alternatively, proceed at your own risk using new_header=create_header(observer,map) and map.reproject_to(new_header) .")
outmap = map.reproject_to(obs_metadict_header)
outmap.plot_settings = map.plot_settings
return outmap

Expand All @@ -212,7 +264,7 @@ def plot_map_reproj(map, reprojected_map):
fig = plt.figure(figsize=(16, 4))
ax1 = fig.add_subplot(1, 3, 1, projection=map)
map.plot(axes=ax1, title=f'Input {map.detector} map {map.date}')
reprojected_map.draw_grid(annotate=False, color='w')
map.draw_grid(annotate=False, color='w')
ax2 = fig.add_subplot(1, 3, 2, projection=reprojected_map)
reprojected_map.plot(axes=ax2, title=f'Map as seen by observer {map.date}')
reprojected_map.draw_grid(annotate=False, color='k')
Expand Down