Overlaying an IRIS SJI "map" with a Sunpy AIA map

[Tiago Pereira]

Hi,

I would like to overlay an IRIS SJI image on top of an AIA image. In the
Sunpy documentation there is an example on how to plot different map
objects, but unfortunately this doesn't work with an
irispy.sji.IRISMapCube object, which is subclassed from ndcube. Example:


>>> my_maps = sunpy.map.Map(aia_map, composite=True)
>>> my_maps.add_map(iris_map[0])
>>> my_maps.set_alpha(1, 0.2)
>>> my_maps.plot()
(...)
AttributeError: 'IRISMapCube' object has no attribute '_get_lon_lat'

It would be nice if something like the above would work one day. I guess
ndcube can have very different types of axes (not just lat/lon), but it
should be possible to do a check.

Anyway, to proceed with my plot I tried to use WCSAxes, but I still got
no success. In my case the IRIS image is rolled, so the image axes do
not match the solar xy. To account for this, I rotated the AIA image
with the inverse WCS rotation matrix of IRIS (visual inspection confirms
the result was fine):

>>> aia_rot = aia_map.rotate(rmatrix=np.matrix(iris_map.wcs.wcs.pc[:-1,
:-1]).I)

Then I start the plot:

>>> ax = plt.subplot(projection=aia_rot.wcs)
>>> ax.imshow(aia_rot.data**0.4)
>>> ax.imshow(iris_map[0].data, cmap='gist_heat', vmin=0,
... transform=ax.get_transform(iris_map.wcs.dropaxis(-1)),
... alpha=0.1)

If I enter the above, it fails with a very long error message (below),
but if I ran either of the imshow commands separately, it works. Here's
the full error message:

---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
~/miniconda/lib/python3.6/site-packages/IPython/core/formatters.py in
__call__(self, obj)
339 pass
340 else:
--> 341 return printer(obj)
342 # Finally look for special method names
343 method = get_real_method(obj, self.print_method)

~/miniconda/lib/python3.6/site-packages/IPython/core/pylabtools.py in
<lambda>(fig)
239
240 if 'png' in formats:
--> 241 png_formatter.for_type(Figure, lambda fig:
print_figure(fig, 'png', **kwargs))
242 if 'retina' in formats or 'png2x' in formats:
243 png_formatter.for_type(Figure, lambda fig:
retina_figure(fig, **kwargs))

~/miniconda/lib/python3.6/site-packages/IPython/core/pylabtools.py in
print_figure(fig, fmt, bbox_inches, **kwargs)
123
124 bytes_io = BytesIO()
--> 125 fig.canvas.print_figure(bytes_io, **kw)
126 data = bytes_io.getvalue()
127 if fmt == 'svg':

~/miniconda/lib/python3.6/site-packages/matplotlib/backend_bases.py in
print_figure(self, filename, dpi, facecolor, edgecolor, orientation,
format, **kwargs)
2210 orientation=orientation,
2211 dryrun=True,
-> 2212 **kwargs)
2213 renderer = self.figure._cachedRenderer
2214 bbox_inches = self.figure.get_tightbbox(renderer)

~/miniconda/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py
in print_png(self, filename_or_obj, *args, **kwargs)
511
512 def print_png(self, filename_or_obj, *args, **kwargs):
--> 513 FigureCanvasAgg.draw(self)
514 renderer = self.get_renderer()
515 original_dpi = renderer.dpi

~/miniconda/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py
in draw(self)
431 # if toolbar:
432 # toolbar.set_cursor(cursors.WAIT)
--> 433 self.figure.draw(self.renderer)
434 # A GUI class may be need to update a window using
this draw, so
435 # don't forget to call the superclass.

~/miniconda/lib/python3.6/site-packages/matplotlib/artist.py 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:

~/miniconda/lib/python3.6/site-packages/matplotlib/figure.py in
draw(self, renderer)
1473
1474 mimage._draw_list_compositing_images(
-> 1475 renderer, self, artists, self.suppressComposite)
1476
1477 renderer.close_group('figure')

~/miniconda/lib/python3.6/site-packages/matplotlib/image.py in
_draw_list_compositing_images(renderer, parent, artists, suppress_composite)
139 if not_composite or not has_images:
140 for a in artists:
--> 141 a.draw(renderer)
142 else:
143 # Composite any adjacent images together

~/miniconda/lib/python3.6/site-packages/astropy/visualization/wcsaxes/core.py
in draw(self, renderer, inframe)
390 self.coords.frame._update_patch_path()
391
--> 392 super().draw(renderer, inframe=inframe)
393
394 self._drawn = True

~/miniconda/lib/python3.6/site-packages/matplotlib/artist.py 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:

~/miniconda/lib/python3.6/site-packages/matplotlib/axes/_base.py in
draw(self, renderer, inframe)
2605 renderer.stop_rasterizing()
2606
-> 2607 mimage._draw_list_compositing_images(renderer, self,
artists)
2608
2609 renderer.close_group('axes')

~/miniconda/lib/python3.6/site-packages/matplotlib/image.py in
_draw_list_compositing_images(renderer, parent, artists, suppress_composite)
139 if not_composite or not has_images:
140 for a in artists:
--> 141 a.draw(renderer)
142 else:
143 # Composite any adjacent images together

~/miniconda/lib/python3.6/site-packages/matplotlib/artist.py 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:

~/miniconda/lib/python3.6/site-packages/matplotlib/image.py in
draw(self, renderer, *args, **kwargs)
591 else:
592 im, l, b, trans = self.make_image(
--> 593 renderer, renderer.get_image_magnification())
594 if im is not None:
595 renderer.draw_image(gc, l, b, im)

~/miniconda/lib/python3.6/site-packages/matplotlib/image.py in
make_image(self, renderer, magnification, unsampled)
839 return self._make_image(
840 self._A, bbox, transformed_bbox, self.axes.bbox,
magnification,
--> 841 unsampled=unsampled)
842
843 def _check_unsampled_image(self, renderer):

~/miniconda/lib/python3.6/site-packages/matplotlib/image.py in
_make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification,
unsampled, round_to_pixel_border)
307 "this method is called.")
308
--> 309 clipped_bbox = Bbox.intersection(out_bbox, clip_bbox)
310
311 if clipped_bbox is None:

~/miniconda/lib/python3.6/site-packages/matplotlib/transforms.py in
intersection(bbox1, bbox2)
757 if they do not intersect.
758 """
--> 759 x0 = np.maximum(bbox1.xmin, bbox2.xmin)
760 x1 = np.minimum(bbox1.xmax, bbox2.xmax)
761 y0 = np.maximum(bbox1.ymin, bbox2.ymin)

~/miniconda/lib/python3.6/site-packages/matplotlib/transforms.py in
xmin(self)
353 :attr:`xmin` is the left edge of the bounding box.
354 """
--> 355 return np.min(self.get_points()[:, 0])
356
357 @property

~/miniconda/lib/python3.6/site-packages/matplotlib/transforms.py in
get_points(self)
1080 [p[1, 0], p[0, 1]],
1081 [p[0, 0], p[1, 1]],
-> 1082 [p[1, 0], p[1, 1]]])
1083 points = np.ma.filled(points, 0.0)
1084

~/miniconda/lib/python3.6/site-packages/matplotlib/transforms.py in
transform(self, values)
1444
1445 # Transform the values
-> 1446 res =
self.transform_affine(self.transform_non_affine(values))
1447
1448 # Convert the result back to the shape of the input values.

~/miniconda/lib/python3.6/site-packages/matplotlib/transforms.py in
transform_non_affine(self, points)
2487 return points
2488 elif not self._a.is_affine and self._b.is_affine:
-> 2489 return self._a.transform_non_affine(points)
2490 else:
2491 return self._b.transform_non_affine(

~/miniconda/lib/python3.6/site-packages/matplotlib/transforms.py in
transform_non_affine(self, points)
2490 else:
2491 return self._b.transform_non_affine(
-> 2492 self._a.transform(points))
2493 transform_non_affine.__doc__ =
Transform.transform_non_affine.__doc__
2494

~/miniconda/lib/python3.6/site-packages/matplotlib/transforms.py in
transform_non_affine(self, points)
2490 else:
2491 return self._b.transform_non_affine(
-> 2492 self._a.transform(points))
2493 transform_non_affine.__doc__ =
Transform.transform_non_affine.__doc__
2494

~/miniconda/lib/python3.6/site-packages/astropy/visualization/wcsaxes/transforms.py
in transform(self, input_coords)
261 # on all values and just ignore Numpy warnings
262 with np.errstate(all='ignore'):
--> 263 c_out = c_in.transform_to(self.output_system)
264
265 if issubclass(c_out.representation,
(SphericalRepresentation, UnitSphericalRepresentation)):

~/miniconda/lib/python3.6/site-packages/astropy/coordinates/sky_coordinate.py
in transform_to(self, frame, merge_attributes)
537 # Do the transformation, returning a coordinate frame of
the desired
538 # final type (not generic).
--> 539 new_coord = trans(self.frame, generic_frame)
540
541 # Finally make the new SkyCoord object from the
`new_coord` and

~/miniconda/lib/python3.6/site-packages/astropy/coordinates/transformations.py
in __call__(self, fromcoord, toframe)
1374
1375 curr_toframe = t.tosys(**frattrs)
-> 1376 curr_coord = t(curr_coord, curr_toframe)
1377
1378 # this is safe even in the case where self.transforms is
empty, because

~/miniconda/lib/python3.6/site-packages/astropy/coordinates/transformations.py
in __call__(self, fromcoord, toframe)
813
814 def __call__(self, fromcoord, toframe):
--> 815 res = self.func(fromcoord, toframe)
816 if not isinstance(res, self.tosys):
817 raise TypeError('the transformation function yielded
{0} but '

~/miniconda/lib/python3.6/site-packages/sunpy/coordinates/transformations.py
in hpc_to_hpc(heliopcoord, heliopframe)
235 from astropy.tests.helper import quantity_allclose
236 if (heliopcoord.observer == heliopframe.observer or
--> 237 (quantity_allclose(heliopcoord.observer.lat,
heliopframe.observer.lat) and
238 quantity_allclose(heliopcoord.observer.lon,
heliopframe.observer.lon) and
239 quantity_allclose(heliopcoord.observer.radius,
heliopframe.observer.radius))):

AttributeError: 'str' object has no attribute 'lat'



Any idea on what is going wrong here?


Thanks,

Tiago