Merge pull request #1352 from Cadair/map_wcs

Add a on-the-fly wcs property for Map

sunpy/map/mapbase.py

@@ -16,6 +16,7 @@
 from matplotlib import patches
 from matplotlib import cm
 
+import astropy.wcs
 from .nddata_compat import NDDataCompat as NDData
 
 from sunpy.image.transform import affine_transform
@@ -163,6 +164,17 @@ def __repr__(self):
            dx=self.scale['x'], dy=self.scale['y'])
 + self.data.__repr__())
 
+    @property
+    def wcs(self):
+        w2 = astropy.wcs.WCS(naxis=2)
+        w2.wcs.crpix = [self.reference_pixel['x'], self.reference_pixel['y']]
+        w2.wcs.cdelt = [self.scale['x'], self.scale['y']]
+        w2.wcs.crval = [self.reference_coordinate['x'], self.reference_coordinate['y']]
+        w2.wcs.ctype = [self.coordinate_system['x'], self.coordinate_system['y']]
+        w2.wcs.pc = self.rotation_matrix
+        w2.wcs.cunit = [self.units['x'], self.units['y']]
+
+        return w2
 
     #Some numpy extraction
     @property

sunpy/map/tests/test_mapbase.py

@@ -10,6 +10,7 @@
 
 import pytest
 
+import astropy.wcs
 from astropy.io import fits
 import astropy.units as u
 
@@ -66,6 +67,21 @@ def test_repr_no_obs(generic_map):
 def test_repr_obs(aia171_test_map):
     assert aia171_test_map.__repr__() == 'SunPy AIAMap\n---------\nObservatory:\t SDO\nInstrument:\t AIA_3\nDetector:\t AIA\nMeasurement:\t 171\nObs Date:\t 2011-02-15T00:00:00.34\ndt:\t\t 2.000191\nDimension:\t [128, 128]\n[dx, dy] =\t [19.183648, 19.183648]\n\narray([[-1.25,  0.  ,  1.  , ...,  0.  ,  0.5 , -0.75],\n       [ 0.75, -0.25, -0.5 , ...,  0.25,  0.  , -0.25],\n       [ 0.  ,  0.5 ,  1.75, ...,  0.  ,  0.5 ,  0.  ],\n       ..., \n       [ 1.  ,  0.25, -0.25, ...,  0.  ,  0.  ,  0.  ],\n       [-0.25,  0.  , -0.5 , ...,  0.75, -0.75,  0.  ],\n       [ 0.75,  1.5 , -0.75, ...,  0.  , -0.5 ,  0.5 ]])'
 
+def test_wcs(aia171_test_map):
+    wcs = aia171_test_map.wcs
+    assert isinstance(wcs, astropy.wcs.WCS)
+
+    assert all(wcs.wcs.crpix == [aia171_test_map.reference_pixel['x'],
+                                 aia171_test_map.reference_pixel['y']])
+    assert all(wcs.wcs.cdelt == [aia171_test_map.scale['x'],
+                                 aia171_test_map.scale['y']])
+    assert all(wcs.wcs.crval == [aia171_test_map.reference_coordinate['x'],
+                                 aia171_test_map.reference_coordinate['y']])
+    assert set(wcs.wcs.ctype) == set([aia171_test_map.coordinate_system['x'],
+                             aia171_test_map.coordinate_system['y']])
+    np.testing.assert_allclose(wcs.wcs.pc, aia171_test_map.rotation_matrix)
+    assert set(wcs.wcs.cunit) == set([u.Unit(a) for a in [aia171_test_map.units['x'],
+                                                          aia171_test_map.units['y']]])
 
 def test_dtype(generic_map):
     assert generic_map.dtype == np.float64
@@ -210,7 +226,7 @@ def test_rotation_matrix_cd_cdelt_square():
 def test_swap_cd():
     amap = sunpy.map.Map(os.path.join(testpath, 'swap_lv1_20140606_000113.fits'))
     np.testing.assert_allclose(amap.rotation_matrix, np.matrix([[1., 0], [0, 1.]]))
-    
+
 
 def test_data_range(generic_map):
     """Make sure xrange and yrange work"""