In [1]:
import sunpy
import os,sys, re
import datetime
import numpy as np
from sunpy.map import Map
from astropy.io import fits as pyfits
from astropy import wcs
from astropy.coordinates import EarthLocation, SkyCoord, Angle
import matplotlib.pyplot as plt
import copy
import astropy.units as u
import datetime
%matplotlib inline


/Users/kkozarev/anaconda2/lib/python2.7/site-packages/glymur/config.py:122: UserWarning: The openjp2 library at /opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/libopenjp2.dylib could not be loaded.
  warnings.warn(msg, UserWarning)

In [18]:
#Create the AIA map
#aiafile='/Users/kkozarev/sunpy/data/sample_data/AIA20110319_105400_0171.fits'
aiafile='/Users/kkozarev/CORWAV_DATA/AIA_data/2015/11/04/H0300/AIA20151104_033847_0171.fits'
aiamap = Map(aiafile)
aiamap.coordinate_frame
aiamap.coordinate_system


Out[18]:
SpatialPair(axis1='HPLN-TAN', axis2='HPLT-TAN')

In [3]:
#Load the MWA data
#infile='/Volumes/PLUME/MWA_DATA/093-094/1130643536/1130643536_c093-094_f8-14_t033843_XX_d001.fits'
infile='/Volumes/PLUME/MWA_DATA/062-063/1130643536/1130643536_c062-063_f8-14_t033843_XX_d001.fits'

fh = pyfits.open(infile)
data=fh[0].data[0,0,:,:]
fh[0].header['ORIGIN'] = 'CASA 4.7.0-REL (r38335)'
print('Header origin: "{0}"'.format(fh[0].header['ORIGIN']))
header = copy.copy(fh[0].header)
header['ORIGIN'] = 'BLANK'
header.update()
try:
    print(wcs.WCS(header))
except ValueError:
    print(wcs.WCS(header))


Header origin: "CASA 4.7.0-REL (r38335)"
WCS Keywords

Number of WCS axes: 4
CTYPE : 'RA---SIN'  'DEC--SIN'  'FREQ'  'STOKES'  
CRVAL : 218.13355931640001  -15.18065681751  79160000.01376  -5.0  
CRPIX : 513.0  513.0  1.0  1.0  
PC1_1 PC1_2 PC1_3 PC1_4  : 1.0  0.0  0.0  0.0  
PC2_1 PC2_2 PC2_3 PC2_4  : 0.0  1.0  0.0  0.0  
PC3_1 PC3_2 PC3_3 PC3_4  : 0.0  0.0  1.0  0.0  
PC4_1 PC4_2 PC4_3 PC4_4  : 0.0  0.0  0.0  1.0  
CDELT : -0.0055555555555559999  0.0055555555555559999  280000.00310550001  1.0  
NAXIS : 1024  1024  1  1
WARNING: FITSFixedWarning: PC01_01 = 1.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC02_01 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC03_01 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC04_01 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC01_02 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC02_02 = 1.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC03_02 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC04_02 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC01_03 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC02_03 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC03_03 = 1.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC04_03 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC01_04 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC02_04 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC03_04 = 0.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC04_04 = 1.000000000000E+00 
indices in parameterized keywords must not have leading zeroes. [astropy.wcs.wcs]

In [50]:
obstime=header['date-obs']
#Let's now get the distance between the Earth and the Sun at the time the photo was taken
dsun = sunpy.coordinates.get_sunearth_distance(obstime)
print dsun

#The size of the Sun in the sky is then
rsun_obs = np.arctan(sunpy.sun.constants.radius / dsun).to('arcsec')
print rsun_obs

#The image plate scale is
plate_scale = [float(header['cdelt1']*3600.),float(header['cdelt2']*3600.)]#*u.arcsec/u.pix
#print plate_scale

#We also need the solar rotation angle
obsloc=(float(header['obsgeo-x']),float(header['obsgeo-y']),float(header['obsgeo-z']))
loc = EarthLocation(x=obsloc[0]*u.m,y=obsloc[1]*u.m,z=obsloc[2]*u.m)
fudge_angle = 0.0 * u.deg # update this in case your camera was not perfectly level.
solar_rotation_angle = sunpy.coordinates.get_sun_orientation(loc, header['date-obs']) + fudge_angle
print solar_rotation_angle

hgln_obs = sunpy.coordinates.get_sun_L0(obstime)
hglt_obs = sunpy.coordinates.get_sun_B0(obstime)
print("{0} {1}".format(hglt_obs, hgln_obs))


#Calculate the starting coords in 'HPLN-TAN'
refra=Angle(header['crval1'],u.deg)
refdec=Angle(header['crval2'],u.deg)
ref_coords = SkyCoord(ra=refra,dec=refdec, obstime=obstime,
              distance=dsun,frame='icrs').transform_to(aiamap.coordinate_frame)

#Create a WCS object and put info in
w = wcs.WCS(naxis=2)
w.wcs.dateobs = header["DATE-OBS"]
w.wcs.crpix = [header['crpix1'],header['crpix2']]
w.wcs.ctype = [header['ctype1'],header['ctype2']]
w.wcs.cunit = ['deg', 'deg']
w.wcs.crval = [header['crval1'],header['crval2']]
w.wcs.cdelt = [header['cdelt1'],header['cdelt2']]
#w.wcs.crval = [ref_coords.Tx.value,ref_coords.Ty.value]
#w.wcs.crval = [0.-ref_coords.Tx.value,0.-ref_coords.Ty.value]
#Assume the image is centered on the Sun!!!
#w.wcs.crval = [0.,0.]
#w.wcs.ctype = ['icrs', 'icrs']

#w.wcs.cdelt =  plate_scale

#Convert to a header which will be used to generate the map
mapheader = dict(w.to_header())
mapheader.update({'CROTA': solar_rotation_angle.to('deg').value})
mapheader.update({'DSUN_OBS': dsun.to('m').value})
#mapheader.update({'HGLN_OBS': hgln_obs.to('deg').value})
#mapheader.update({'HGLT_OBS': hglt_obs.to('deg').value})
#mapheader.update({'CTYPE1': 'RA---SIN'})
#mapheader.update({'CTYPE2': 'DEC--SIN'})
mapheader.update({'RSUN': dsun.to('m').value})
mapheader.update({'RSUN_OBS': np.arctan(sunpy.sun.constants.radius / dsun).to('arcsec').value})

mapheader.update({'AUTHOR': header['observer']})
#mapheader.update({'EXPTIME': exposure_time.to('s').value})
mapheader.update({'TELESCOP': header['telescop']})

#Create the MWA map
mwamap = Map(data, mapheader)


0.991890104238 AU
966.796596859 arcsec
175d33m30.4524s
4.1076905811 deg 316.654003591 deg

In [37]:
mwamap


/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:632: Warning: Missing metadata for heliographic latitude: assuming Earth-based observer
  [ 1.    ,  0.1875, -0.8125,  0.125 ,  0.3125],
/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:632: Warning: Missing metadata for heliographic latitude: assuming Earth-based observer
  [ 1.    ,  0.1875, -0.8125,  0.125 ,  0.3125],
/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:632: Warning: Missing metadata for heliographic latitude: assuming Earth-based observer
  [ 1.    ,  0.1875, -0.8125,  0.125 ,  0.3125],
Out[37]:
SunPy Map
---------
Observatory:		 MWA
Instrument:		 
Detector:		 
Measurement:		 0.0
Wavelength:		 0.0
Observation Date:	 2015-11-04 03:38:41
Exposure Time:		 0.000000 s
Dimension:		 [ 1024.  1024.] pix
Coordinate System:	 icrs
Scale:			 [-0.00555556  0.00555556] deg / pix
Reference Pixel:	 [ 513.  513.] pix
Reference Coord:	 [ 218.13355932  -15.18065682] deg

array([[-0.3892709 , -0.41520908, -0.41844285, ...,  0.27955323,
         0.19957899,  0.11942415],
       [-0.39164123, -0.4053531 , -0.39761546, ...,  0.26869285,
         0.19014969,  0.11316711],
       [-0.39393985, -0.3958348 , -0.37788472, ...,  0.24971271,
         0.17237014,  0.09819561],
       ..., 
       [-0.69270629, -0.72989023, -0.73872876, ...,  0.1622441 ,
         0.16237251,  0.16370109],
       [-0.72364223, -0.75563061, -0.75918961, ...,  0.22050744,
         0.22246408,  0.22616953],
       [-0.75068754, -0.77595532, -0.77301866, ...,  0.26766422,
         0.26971716,  0.27415675]], dtype=float32)

In [163]:
fig = plt.figure(figsize=(10,10))
ax = plt.subplot(projection=aiamap)
aiamap.plot(axes=ax)
aiamap.draw_grid(axes=ax)
aiamap.draw_limb(axes=ax)
plt.show()



In [51]:
fig = plt.figure(figsize=(10,10))
ax = plt.subplot(projection=mwamap)
mwamap.plot(axes=ax)
mwamap.draw_grid(axes=ax)
mwamap.draw_limb(axes=ax)
plt.show()


/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:632: Warning: Missing metadata for heliographic latitude: assuming Earth-based observer
  [ 1.    ,  0.1875, -0.8125,  0.125 ,  0.3125],
---------------------------------------------------------------------------
UnitsError                                Traceback (most recent call last)
/Users/kkozarev/anaconda2/lib/python2.7/site-packages/IPython/core/formatters.pyc in __call__(self, obj)
    305                 pass
    306             else:
--> 307                 return printer(obj)
    308             # Finally look for special method names
    309             method = get_real_method(obj, self.print_method)

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in <lambda>(fig)
    225 
    226     if 'png' in formats:
--> 227         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    228     if 'retina' in formats or 'png2x' in formats:
    229         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in print_figure(fig, fmt, bbox_inches, **kwargs)
    117 
    118     bytes_io = BytesIO()
--> 119     fig.canvas.print_figure(bytes_io, **kw)
    120     data = bytes_io.getvalue()
    121     if fmt == 'svg':

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2206                     orientation=orientation,
   2207                     dryrun=True,
-> 2208                     **kwargs)
   2209                 renderer = self.figure._cachedRenderer
   2210                 bbox_inches = self.figure.get_tightbbox(renderer)

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
    505 
    506     def print_png(self, filename_or_obj, *args, **kwargs):
--> 507         FigureCanvasAgg.draw(self)
    508         renderer = self.get_renderer()
    509         original_dpi = renderer.dpi

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in draw(self)
    428             if toolbar:
    429                 toolbar.set_cursor(cursors.WAIT)
--> 430             self.figure.draw(self.renderer)
    431         finally:
    432             if toolbar:

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     53                 renderer.start_filter()
     54 
---> 55             return draw(artist, renderer, *args, **kwargs)
     56         finally:
     57             if artist.get_agg_filter() is not None:

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/figure.pyc in draw(self, renderer)
   1293 
   1294             mimage._draw_list_compositing_images(
-> 1295                 renderer, self, artists, self.suppressComposite)
   1296 
   1297             renderer.close_group('figure')

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/image.pyc in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    136     if not_composite or not has_images:
    137         for a in artists:
--> 138             a.draw(renderer)
    139     else:
    140         # Composite any adjacent images together

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/core.pyc in draw(self, renderer, inframe)
    338             coords.frame.update()
    339             for coord in coords:
--> 340                 coord._draw_grid(renderer)
    341 
    342         for coords in self._all_coords:

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/coordinate_helpers.pyc in _draw_grid(self, renderer)
    432         renderer.open_group('grid lines')
    433 
--> 434         self._update_ticks()
    435 
    436         if self.grid_lines_kwargs['visible']:

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/coordinate_helpers.pyc in _update_ticks(self)
    487 
    488         # Find the range of coordinates in all directions
--> 489         coord_range = self.parent_map.get_coord_range()
    490 
    491         # First find the ticks we want to show

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/coordinates_map.pyc in get_coord_range(self)
    161                                      [xmin, xmax, ymin, ymax],
    162                                      [coord.coord_type for coord in self],
--> 163                                      [coord.coord_unit for coord in self])

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/coordinate_range.pyc in find_coordinate_range(transform, extent, coord_types, coord_units)
     45     y = np.linspace(extent[2], extent[3], ny + 1)
     46     xp, yp = np.meshgrid(x, y)
---> 47     world = transform.transform(np.vstack([xp.ravel(), yp.ravel()]).transpose())
     48 
     49     ranges = []

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/transforms.pyc in transform(self, values)
   1421 
   1422         # Transform the values
-> 1423         res = self.transform_affine(self.transform_non_affine(values))
   1424 
   1425         # Convert the result back to the shape of the input values.

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/transforms.pyc in transform_non_affine(self, points)
   2458         else:
   2459             return self._b.transform_non_affine(
-> 2460                                 self._a.transform(points))
   2461     transform_non_affine.__doc__ = Transform.transform_non_affine.__doc__
   2462 

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/transforms.pyc in transform(self, input_coords)
    249         # on all values and just ignore Numpy warnings
    250         with np.errstate(all='ignore'):
--> 251             c_out = c_in.transform_to(self.output_system)
    252 
    253         if issubclass(c_out.representation, (SphericalRepresentation, UnitSphericalRepresentation)):

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/coordinates/sky_coordinate.pyc in transform_to(self, frame, merge_attributes)
    479         # Do the transformation, returning a coordinate frame of the desired
    480         # final type (not generic).
--> 481         new_coord = trans(self.frame, generic_frame)
    482 
    483         # Finally make the new SkyCoord object from the `new_coord` and

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/coordinates/transformations.pyc in __call__(self, fromcoord, toframe)
   1312 
   1313             curr_toframe = t.tosys(**frattrs)
-> 1314             curr_coord = t(curr_coord, curr_toframe)
   1315 
   1316         # this is safe even in the case where self.transforms is empty, because

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/coordinates/transformations.pyc in __call__(self, fromcoord, toframe)
   1127     def __call__(self, fromcoord, toframe):
   1128 
-> 1129         M, vec = self.transform_func(fromcoord, toframe)
   1130         newrep = self._apply_transform(fromcoord, M, vec)
   1131 

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/coordinates/builtin_frames/icrs_cirs_transforms.pyc in icrs_to_hcrs(icrs_coo, hcrs_frame)
    310     # this is just an origin translation so without a distance it cannot go ahead
    311     if isinstance(icrs_coo.data, UnitSphericalRepresentation):
--> 312         raise u.UnitsError(_NEED_ORIGIN_HINT.format(icrs_coo.__class__.__name__))
    313 
    314     if icrs_coo.data.differentials:

UnitsError: The input ICRS coordinates do not have length units. This probably means you created coordinates with lat/lon but no distance.  Heliocentric<->ICRS transforms cannot function in this case because there is an origin shift.
<matplotlib.figure.Figure at 0x114948d90>

In [30]:
#Supposedly de-rotated map
mwamap2 = mwamap.rotate(rmatrix=np.linalg.inv(aiamap.rotation_matrix),
                        recenter=True, order=3, scale=(mwamap.scale[0]/aiamap.scale[0]))#/10


/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:632: Warning: Missing metadata for heliographic latitude: assuming Earth-based observer
  [ 1.    ,  0.1875, -0.8125,  0.125 ,  0.3125],
---------------------------------------------------------------------------
InvalidTransformError                     Traceback (most recent call last)
<ipython-input-30-73ef25806b74> in <module>()
      1 #Supposedly de-rotated map
      2 mwamap2 = mwamap.rotate(rmatrix=np.linalg.inv(aiamap.rotation_matrix),
----> 3                         recenter=True, order=3, scale=(mwamap.scale[0]/aiamap.scale[0]))#/10

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc in rotate(self, angle, rmatrix, order, scale, recenter, missing, use_scipy)
   1060         # The FITS-WCS transform is by definition defined around the
   1061         # reference coordinate in the header.
-> 1062         lon, lat = self._get_lon_lat(self.reference_coordinate.frame)
   1063         rotation_center = u.Quantity([lon, lat])
   1064 

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc in reference_coordinate(self)
    654         return SkyCoord(self._reference_longitude,
    655                         self._reference_latitude,
--> 656                         frame=self.coordinate_frame)
    657 
    658     @property

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc in coordinate_frame(self)
    295         information for this Map.
    296         """
--> 297         return astropy.wcs.utils.wcs_to_celestial_frame(self.wcs)
    298 
    299     def _as_mpl_axes(self):

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/wcs/utils.pyc in wcs_to_celestial_frame(wcs)
    143     for mapping_set in WCS_FRAME_MAPPINGS:
    144         for func in mapping_set:
--> 145             frame = func(wcs)
    146             if frame is not None:
    147                 return frame

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/wcs/utils.pyc in _wcs_to_celestial_frame_builtin(wcs)
     49 
     50     # Keep only the celestial part of the axes
---> 51     wcs = wcs.sub([WCSSUB_CELESTIAL])
     52 
     53     if wcs.wcs.lng == -1 or wcs.wcs.lat == -1:

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/wcs/wcs.pyc in sub(self, axes)
    553 
    554     def sub(self, axes=None):
--> 555         copy = self.deepcopy()
    556         copy.wcs = self.wcs.sub(axes)
    557         copy.naxis = copy.wcs.naxis

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/wcs/wcs.pyc in deepcopy(self)
    550         :mod:`copy` stdlib module.
    551         """
--> 552         return copy.deepcopy(self)
    553 
    554     def sub(self, axes=None):

/Users/kkozarev/anaconda2/lib/python2.7/copy.pyc in deepcopy(x, memo, _nil)
    172             copier = getattr(x, "__deepcopy__", None)
    173             if copier:
--> 174                 y = copier(memo)
    175             else:
    176                 reductor = dispatch_table.get(cls)

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/wcs/wcs.pyc in __deepcopy__(self, memo)
    523                          (deepcopy(self.cpdis1, memo),
    524                           deepcopy(self.cpdis2, memo)),
--> 525                          deepcopy(self.wcs, memo),
    526                          (deepcopy(self.det2im1, memo),
    527                           deepcopy(self.det2im2, memo)))

/Users/kkozarev/anaconda2/lib/python2.7/copy.pyc in deepcopy(x, memo, _nil)
    172             copier = getattr(x, "__deepcopy__", None)
    173             if copier:
--> 174                 y = copier(memo)
    175             else:
    176                 reductor = dispatch_table.get(cls)

InvalidTransformError: ERROR 7 in wcsset() at line 2101 of file cextern/wcslib/C/wcs.c:
Ill-conditioned coordinate transformation parameter.
ERROR 4 in celset() at line 422 of file cextern/wcslib/C/cel.c:
Ill-conditioned coordinate transformation parameters
No valid solution for latp for these values of phip, phi0, and theta0.

In [25]:
#Plot the de-rotated map
fig = plt.figure(figsize=(10,10))
mwamap2.plot()
mwamap2.draw_grid()
mwamap2.draw_limb()
plt.show()


/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:632: Warning: Missing metadata for heliographic latitude: assuming Earth-based observer
  [ 1.    ,  0.1875, -0.8125,  0.125 ,  0.3125],
/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:632: Warning: Missing metadata for heliographic latitude: assuming Earth-based observer
  [ 1.    ,  0.1875, -0.8125,  0.125 ,  0.3125],
/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:632: Warning: Missing metadata for heliographic latitude: assuming Earth-based observer
  [ 1.    ,  0.1875, -0.8125,  0.125 ,  0.3125],
---------------------------------------------------------------------------
UnitsError                                Traceback (most recent call last)
/Users/kkozarev/anaconda2/lib/python2.7/site-packages/IPython/core/formatters.pyc in __call__(self, obj)
    305                 pass
    306             else:
--> 307                 return printer(obj)
    308             # Finally look for special method names
    309             method = get_real_method(obj, self.print_method)

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in <lambda>(fig)
    225 
    226     if 'png' in formats:
--> 227         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    228     if 'retina' in formats or 'png2x' in formats:
    229         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in print_figure(fig, fmt, bbox_inches, **kwargs)
    117 
    118     bytes_io = BytesIO()
--> 119     fig.canvas.print_figure(bytes_io, **kw)
    120     data = bytes_io.getvalue()
    121     if fmt == 'svg':

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2206                     orientation=orientation,
   2207                     dryrun=True,
-> 2208                     **kwargs)
   2209                 renderer = self.figure._cachedRenderer
   2210                 bbox_inches = self.figure.get_tightbbox(renderer)

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
    505 
    506     def print_png(self, filename_or_obj, *args, **kwargs):
--> 507         FigureCanvasAgg.draw(self)
    508         renderer = self.get_renderer()
    509         original_dpi = renderer.dpi

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in draw(self)
    428             if toolbar:
    429                 toolbar.set_cursor(cursors.WAIT)
--> 430             self.figure.draw(self.renderer)
    431         finally:
    432             if toolbar:

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     53                 renderer.start_filter()
     54 
---> 55             return draw(artist, renderer, *args, **kwargs)
     56         finally:
     57             if artist.get_agg_filter() is not None:

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/figure.pyc in draw(self, renderer)
   1293 
   1294             mimage._draw_list_compositing_images(
-> 1295                 renderer, self, artists, self.suppressComposite)
   1296 
   1297             renderer.close_group('figure')

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/image.pyc in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    136     if not_composite or not has_images:
    137         for a in artists:
--> 138             a.draw(renderer)
    139     else:
    140         # Composite any adjacent images together

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/core.pyc in draw(self, renderer, inframe)
    338             coords.frame.update()
    339             for coord in coords:
--> 340                 coord._draw_grid(renderer)
    341 
    342         for coords in self._all_coords:

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/coordinate_helpers.pyc in _draw_grid(self, renderer)
    432         renderer.open_group('grid lines')
    433 
--> 434         self._update_ticks()
    435 
    436         if self.grid_lines_kwargs['visible']:

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/coordinate_helpers.pyc in _update_ticks(self)
    487 
    488         # Find the range of coordinates in all directions
--> 489         coord_range = self.parent_map.get_coord_range()
    490 
    491         # First find the ticks we want to show

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/coordinates_map.pyc in get_coord_range(self)
    161                                      [xmin, xmax, ymin, ymax],
    162                                      [coord.coord_type for coord in self],
--> 163                                      [coord.coord_unit for coord in self])

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/coordinate_range.pyc in find_coordinate_range(transform, extent, coord_types, coord_units)
     45     y = np.linspace(extent[2], extent[3], ny + 1)
     46     xp, yp = np.meshgrid(x, y)
---> 47     world = transform.transform(np.vstack([xp.ravel(), yp.ravel()]).transpose())
     48 
     49     ranges = []

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/transforms.pyc in transform(self, values)
   1421 
   1422         # Transform the values
-> 1423         res = self.transform_affine(self.transform_non_affine(values))
   1424 
   1425         # Convert the result back to the shape of the input values.

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/matplotlib/transforms.pyc in transform_non_affine(self, points)
   2458         else:
   2459             return self._b.transform_non_affine(
-> 2460                                 self._a.transform(points))
   2461     transform_non_affine.__doc__ = Transform.transform_non_affine.__doc__
   2462 

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/visualization/wcsaxes/transforms.pyc in transform(self, input_coords)
    249         # on all values and just ignore Numpy warnings
    250         with np.errstate(all='ignore'):
--> 251             c_out = c_in.transform_to(self.output_system)
    252 
    253         if issubclass(c_out.representation, (SphericalRepresentation, UnitSphericalRepresentation)):

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/coordinates/sky_coordinate.pyc in transform_to(self, frame, merge_attributes)
    479         # Do the transformation, returning a coordinate frame of the desired
    480         # final type (not generic).
--> 481         new_coord = trans(self.frame, generic_frame)
    482 
    483         # Finally make the new SkyCoord object from the `new_coord` and

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/coordinates/transformations.pyc in __call__(self, fromcoord, toframe)
   1312 
   1313             curr_toframe = t.tosys(**frattrs)
-> 1314             curr_coord = t(curr_coord, curr_toframe)
   1315 
   1316         # this is safe even in the case where self.transforms is empty, because

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/coordinates/transformations.pyc in __call__(self, fromcoord, toframe)
   1127     def __call__(self, fromcoord, toframe):
   1128 
-> 1129         M, vec = self.transform_func(fromcoord, toframe)
   1130         newrep = self._apply_transform(fromcoord, M, vec)
   1131 

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/astropy/coordinates/builtin_frames/icrs_cirs_transforms.pyc in icrs_to_hcrs(icrs_coo, hcrs_frame)
    310     # this is just an origin translation so without a distance it cannot go ahead
    311     if isinstance(icrs_coo.data, UnitSphericalRepresentation):
--> 312         raise u.UnitsError(_NEED_ORIGIN_HINT.format(icrs_coo.__class__.__name__))
    313 
    314     if icrs_coo.data.differentials:

UnitsError: The input ICRS coordinates do not have length units. This probably means you created coordinates with lat/lon but no distance.  Heliocentric<->ICRS transforms cannot function in this case because there is an origin shift.
<matplotlib.figure.Figure at 0x118791850>

In [219]:
#Create a composite map for overlaying the radio observations on AIA image
mymaps = Map(aiamap,composite=True)
#Add the de-rotated map
mymaps.add_map(mwamap2)
mymaps.set_levels(1, [50, 60, 70, 80, 90,99], percent = True)

fig = plt.figure(figsize=(10,10))
mymaps.plot()
plt.show()



In [155]:
sunpy.coordinates.get_sun_orientation(loc, obstime)


Out[155]:
$175^\circ33{}^\prime30.4524{}^{\prime\prime}$

In [ ]:
#def remove_non_ascii(text):
#    return unidecode(unicode(text, encoding = "utf-8"))

def mwa_prep(infile):
    data, header = pyfits.getdata(infile, header=True)
    data=data[0,0,:,:]
    w = wcs.WCS(naxis=2)
    w.wcs.naxis1=header["NAXIS1"]
    w.wcs.naxis2=header["NAXIS2"]
    w.wcs.dateobs=header["DATE-OBS"]
    w.wcs.ctype1=header['ctype1']
    w.wcs.ctype2=header['ctype2']
    w.wcs.crpix1=header['crpix1']
    w.wcs.crpix2=header['crpix2']
    w.wcs.crval1=header['crval1']
    w.wcs.crval2=header['crval2']
    w.wcs.cdelt=[header['cdelt1'],header['cdelt2']]