Repository: sunpy/sunpy

Issue:
`.resample()` and `.superpixel()` modify `CDELTi`/`PCi_j` in a way that results in non-intuitive behavior
The FITS WCS standard allows for "considerable flexibility" in the way that `CDELTi` values and the `PCi_j` matrix are defined.  You can choose any combination that multiplies to the `CDi_j` matrix that you need.  However, for 2D imagery, most users would probably expect that the `PCi_j` matrix be as close to a pure rotation matrix as possible, and the `CDELTi` values to be represent the pixel plate scales.

For whatever reason, when modifying a map that uses the `PCi_j` matrix, both `.resample()` and `.superpixel()` modify only the `PCi_j` matrix and do not modify the `CDELTi` values.  Again, this is technically permitted.  However, this has several non-intuitive consequences:

- `.scale` returns `CDELTi` values, which are no longer the pixel plate scales
- `.rotation_matrix` returns the `PCi_j` matrix, which is no longer close to a pure rotation matrix and contains scaling factors
- A de-rotating `.rotate()` call on the map leads to an unexpected resampling because it uses `.rotation_matrix`, which now has scaling factors in it

I'll work up a PR to fix this.

Discussion:

Patch:
diff --git a/changelog/6571.bugfix.rst b/changelog/6571.bugfix.rst
new file mode 100644
index 000000000..3a205e442
--- /dev/null
+++ b/changelog/6571.bugfix.rst
@@ -0,0 +1,2 @@
+Fixed the metadata handling of :meth:`~sunpy.map.GenericMap.resample` and :meth:`~sunpy.map.GenericMap.superpixel` so that the CDELTi values are scaled and the PCi_j matrix (if used) is modified in the correct manner for asymmetric scaling.
+The previous approach of having the PCi_j matrix store all of the scaling resulted in non-intuitive behaviors when accessing the `~sunpy.map.GenericMap.scale` and `~sunpy.map.GenericMap.rotation_matrix` properties, and when de-rotating a map via :meth:`~sunpy.map.GenericMap.rotate`.
diff --git a/sunpy/map/mapbase.py b/sunpy/map/mapbase.py
index 26a960521..1eb448816 100644
--- a/sunpy/map/mapbase.py
+++ b/sunpy/map/mapbase.py
@@ -1542,20 +1542,14 @@ class GenericMap(NDData):
         new_meta = self.meta.copy()
 
         # Update metadata
+        for key in {'cdelt1', 'cd1_1', 'cd2_1'} & self.meta.keys():
+            new_meta[key] *= scale_factor_x
+        for key in {'cdelt2', 'cd1_2', 'cd2_2'} & self.meta.keys():
+            new_meta[key] *= scale_factor_y
         if 'pc1_1' in self.meta:
-            new_meta['pc1_1'] *= scale_factor_x
-            new_meta['pc2_1'] *= scale_factor_x
-            new_meta['pc1_2'] *= scale_factor_y
-            new_meta['pc2_2'] *= scale_factor_y
-        if 'cd1_1' in self.meta:
-            new_meta['cd1_1'] *= scale_factor_x
-            new_meta['cd2_1'] *= scale_factor_x
-            new_meta['cd1_2'] *= scale_factor_y
-            new_meta['cd2_2'] *= scale_factor_y
-        if 'cd1_1' not in self.meta and 'pc1_1' not in self.meta:
-            # Using the CROTA2 and CDELT formalism
-            new_meta['cdelt1'] *= scale_factor_x
-            new_meta['cdelt2'] *= scale_factor_y
+            new_meta['pc1_2'] *= scale_factor_y / scale_factor_x
+            new_meta['pc2_1'] *= scale_factor_x / scale_factor_y
+
         new_meta['crpix1'] = (self.reference_pixel.x.to_value(u.pix) + 0.5) / scale_factor_x + 0.5
         new_meta['crpix2'] = (self.reference_pixel.y.to_value(u.pix) + 0.5) / scale_factor_y + 0.5
         new_meta['naxis1'] = new_data.shape[1]
@@ -2052,25 +2046,20 @@ class GenericMap(NDData):
         new_meta = self.meta.copy()
 
         # Update metadata
+        for key in {'cdelt1', 'cd1_1', 'cd2_1'} & self.meta.keys():
+            new_meta[key] *= dimensions[0]
+        for key in {'cdelt2', 'cd1_2', 'cd2_2'} & self.meta.keys():
+            new_meta[key] *= dimensions[1]
         if 'pc1_1' in self.meta:
-            new_meta['pc1_1'] *= dimensions[0]
-            new_meta['pc2_1'] *= dimensions[0]
-            new_meta['pc1_2'] *= dimensions[1]
-            new_meta['pc2_2'] *= dimensions[1]
-        if 'cd1_1' in self.meta:
-            new_meta['cd1_1'] *= dimensions[0]
-            new_meta['cd2_1'] *= dimensions[0]
-            new_meta['cd1_2'] *= dimensions[1]
-            new_meta['cd2_2'] *= dimensions[1]
-        if 'cd1_1' not in self.meta and 'pc1_1' not in self.meta:
-            # Using the CROTA2 and CDELT formalism
-            new_meta['cdelt1'] *= dimensions[0]
-            new_meta['cdelt2'] *= dimensions[1]
+            new_meta['pc1_2'] *= dimensions[1] / dimensions[0]
+            new_meta['pc2_1'] *= dimensions[0] / dimensions[1]
 
         new_meta['crpix1'] = ((self.reference_pixel.x.to_value(u.pix) +
                                0.5 - offset[0]) / dimensions[0]) + 0.5
         new_meta['crpix2'] = ((self.reference_pixel.y.to_value(u.pix) +
                                0.5 - offset[1]) / dimensions[1]) + 0.5
+        new_meta['naxis1'] = new_array.shape[1]
+        new_meta['naxis2'] = new_array.shape[0]
 
         # Create new map instance
         if self.mask is not None:


Test patch:
diff --git a/sunpy/map/tests/test_mapbase.py b/sunpy/map/tests/test_mapbase.py
index 0cb96b580..d83e83b13 100644
--- a/sunpy/map/tests/test_mapbase.py
+++ b/sunpy/map/tests/test_mapbase.py
@@ -686,8 +686,17 @@ def test_resample_metadata(generic_map, sample_method, new_dimensions):
     Check that the resampled map has correctly adjusted metadata.
     """
     resampled_map = generic_map.resample(new_dimensions, method=sample_method)
-    assert resampled_map.meta['cdelt1'] == generic_map.meta['cdelt1']
-    assert resampled_map.meta['cdelt2'] == generic_map.meta['cdelt2']
+
+    scale_x = generic_map.data.shape[1] / resampled_map.data.shape[1]
+    scale_y = generic_map.data.shape[0] / resampled_map.data.shape[0]
+
+    assert resampled_map.meta['cdelt1'] == scale_x * generic_map.meta['cdelt1']
+    assert resampled_map.meta['cdelt2'] == scale_y * generic_map.meta['cdelt2']
+    assert resampled_map.meta['pc1_1'] == generic_map.meta['pc1_1']
+    assert resampled_map.meta['pc1_2'] == scale_y / scale_x * generic_map.meta['pc1_2']
+    assert resampled_map.meta['pc2_1'] == scale_x / scale_y * generic_map.meta['pc2_1']
+    assert resampled_map.meta['pc2_2'] == generic_map.meta['pc2_2']
+
     # TODO: we should really test the numbers here, not just that the correct
     # header values have been modified. However, I am lazy and we have figure
     # tests.
@@ -698,8 +707,8 @@ def test_resample_metadata(generic_map, sample_method, new_dimensions):
     assert resampled_map.meta['naxis1'] == new_dimensions[0].value
     assert resampled_map.meta['naxis2'] == new_dimensions[1].value
     for key in generic_map.meta:
-        if key not in ('crpix1', 'crpix2', 'crval1',
-                       'crval2', 'naxis1', 'naxis2') and not key.startswith('pc'):
+        if key not in ('cdelt1', 'cdelt2', 'pc1_2', 'pc2_1', 'crpix1', 'crpix2', 'crval1',
+                       'crval2', 'naxis1', 'naxis2'):
             assert resampled_map.meta[key] == generic_map.meta[key]
 
 
@@ -753,6 +762,32 @@ def test_superpixel_dims_values(aia171_test_map, f):
     assert_quantity_allclose(superpix_map.data[0, 0], expected)
 
 
+@pytest.mark.parametrize('f, dimensions', [(np.sum, (2, 3)*u.pix),
+                                           (np.mean, (3, 2)*u.pix)])
+def test_superpixel_metadata(generic_map, f, dimensions):
+    superpix_map = generic_map.superpixel(dimensions, func=f)
+
+    scale_x, scale_y = dimensions.value
+
+    assert superpix_map.meta['cdelt1'] == scale_x * generic_map.meta['cdelt1']
+    assert superpix_map.meta['cdelt2'] == scale_y * generic_map.meta['cdelt2']
+    assert superpix_map.meta['pc1_1'] == generic_map.meta['pc1_1']
+    assert superpix_map.meta['pc1_2'] == scale_y / scale_x * generic_map.meta['pc1_2']
+    assert superpix_map.meta['pc2_1'] == scale_x / scale_y * generic_map.meta['pc2_1']
+    assert superpix_map.meta['pc2_2'] == generic_map.meta['pc2_2']
+
+    assert superpix_map.meta['crpix1'] - 0.5 == (generic_map.meta['crpix1'] - 0.5) / scale_x
+    assert superpix_map.meta['crpix2'] - 0.5 == (generic_map.meta['crpix2'] - 0.5) / scale_y
+    assert u.allclose(superpix_map.meta['crval1'], generic_map.meta['crval1'])
+    assert u.allclose(superpix_map.meta['crval2'], generic_map.meta['crval2'])
+    assert superpix_map.meta['naxis1'] == generic_map.meta['naxis1'] / scale_x
+    assert superpix_map.meta['naxis2'] == generic_map.meta['naxis2'] / scale_y
+    for key in generic_map.meta:
+        if key not in ('cdelt1', 'cdelt2', 'pc1_2', 'pc2_1', 'crpix1', 'crpix2', 'crval1',
+                       'crval2', 'naxis1', 'naxis2'):
+            assert superpix_map.meta[key] == generic_map.meta[key]
+
+
 def test_superpixel_masked(aia171_test_map_with_mask):
     input_dims = u.Quantity(aia171_test_map_with_mask.dimensions)
     dimensions = (2, 2) * u.pix
@@ -817,6 +852,8 @@ def test_resample_rotated_map_cd(cd, method):
     smap = make_simple_map()
 
     smap.meta.update(cd)
+    for key in ['cdelt1', 'cdelt2', 'pc1_1', 'pc1_2', 'pc2_1', 'pc2_2']:
+        del smap.meta[key]
     # Check superpixel with a rotated map with unequal resampling
     new_dims = (1, 2) * u.pix
     new_map = getattr(smap, method)(new_dims)
@@ -1428,6 +1465,7 @@ def test_derotating_nonpurerotation_pcij(aia171_test_map, method):
 
     ax2 = fig.add_subplot(122, projection=derotated_map)
     derotated_map.plot(axes=ax2, title=f'De-rotated map via {method}')
+    ax2.set_aspect(derotated_map.scale[1] / derotated_map.scale[0])
 
     return fig
 
diff --git a/sunpy/tests/figure_hashes_mpl_352_ft_261_astropy_51_animators_100.json b/sunpy/tests/figure_hashes_mpl_352_ft_261_astropy_51_animators_100.json
index a487b30bb..7ddbd55f8 100644
--- a/sunpy/tests/figure_hashes_mpl_352_ft_261_astropy_51_animators_100.json
+++ b/sunpy/tests/figure_hashes_mpl_352_ft_261_astropy_51_animators_100.json
@@ -9,9 +9,9 @@
   "sunpy.map.tests.test_compositemap.test_set_alpha_composite_map": "f811cf6be8abfcd1b3a7f4bae7781806b37b1b2b0cb77ac34c4ed4491a1d633a",
   "sunpy.map.tests.test_compositemap.test_peek_composite_map": "d459182443472330df9b0dde3bde647479c100f56f45b6321619fbf67810c174",
   "sunpy.map.tests.test_mapbase.test_rotation_rect_pixelated_data": "5b37819c36a50a2a7a26944a2f1481027dd17390e73158f6b5510034102d35bb",
-  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[scipy]": "4184d5536f49d2811803ba831e212a7708046f4cb2c9fca1252a017fa9d90127",
-  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[scikit-image]": "27862d38a9ab0c5e61229737e5fe0af962d94ee421e7d5acd2e14c7018bdde44",
-  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[opencv]": "369c75feab7a0b1ebe65bd0ec9084b3c941db45445d88c667311d81f01101afe",
+  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[scipy]": "9168fca9a4e04def3da12624196b29ebe954004ad269faa4bab11276a85366af",
+  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[scikit-image]": "17748c119a97c112e77fcb368af595f2b5faada067a7188abe846c7b160f9ea8",
+  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[opencv]": "cd86fb9c0ed0e26eaae00d816d0b5d3f381b2be22b8e99d45afe5a511ccd55e5",
   "sunpy.map.tests.test_mapsequence.test_norm_animator": "5a7217dd55ab3a593a27741fecfe0851208fd1d31baeca4cbd018d65aaf3f483",
   "sunpy.map.tests.test_mapsequence.test_map_sequence_plot": "3cdc53f41914028a4e9f23112d7e616ce8a423ef5670c8365fdadbc0e59780ca",
   "sunpy.map.tests.test_mapsequence.test_map_sequence_plot_custom_cmap_norm": "c1276702147f4debaa14c17f6b0b2ab7f9e0b1ace309c72e767b2af071c4405e",
diff --git a/sunpy/tests/figure_hashes_mpl_dev_ft_261_astropy_dev_animators_dev.json b/sunpy/tests/figure_hashes_mpl_dev_ft_261_astropy_dev_animators_dev.json
index a487b30bb..7ddbd55f8 100644
--- a/sunpy/tests/figure_hashes_mpl_dev_ft_261_astropy_dev_animators_dev.json
+++ b/sunpy/tests/figure_hashes_mpl_dev_ft_261_astropy_dev_animators_dev.json
@@ -9,9 +9,9 @@
   "sunpy.map.tests.test_compositemap.test_set_alpha_composite_map": "f811cf6be8abfcd1b3a7f4bae7781806b37b1b2b0cb77ac34c4ed4491a1d633a",
   "sunpy.map.tests.test_compositemap.test_peek_composite_map": "d459182443472330df9b0dde3bde647479c100f56f45b6321619fbf67810c174",
   "sunpy.map.tests.test_mapbase.test_rotation_rect_pixelated_data": "5b37819c36a50a2a7a26944a2f1481027dd17390e73158f6b5510034102d35bb",
-  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[scipy]": "4184d5536f49d2811803ba831e212a7708046f4cb2c9fca1252a017fa9d90127",
-  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[scikit-image]": "27862d38a9ab0c5e61229737e5fe0af962d94ee421e7d5acd2e14c7018bdde44",
-  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[opencv]": "369c75feab7a0b1ebe65bd0ec9084b3c941db45445d88c667311d81f01101afe",
+  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[scipy]": "9168fca9a4e04def3da12624196b29ebe954004ad269faa4bab11276a85366af",
+  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[scikit-image]": "17748c119a97c112e77fcb368af595f2b5faada067a7188abe846c7b160f9ea8",
+  "sunpy.map.tests.test_mapbase.test_derotating_nonpurerotation_pcij[opencv]": "cd86fb9c0ed0e26eaae00d816d0b5d3f381b2be22b8e99d45afe5a511ccd55e5",
   "sunpy.map.tests.test_mapsequence.test_norm_animator": "5a7217dd55ab3a593a27741fecfe0851208fd1d31baeca4cbd018d65aaf3f483",
   "sunpy.map.tests.test_mapsequence.test_map_sequence_plot": "3cdc53f41914028a4e9f23112d7e616ce8a423ef5670c8365fdadbc0e59780ca",
   "sunpy.map.tests.test_mapsequence.test_map_sequence_plot_custom_cmap_norm": "c1276702147f4debaa14c17f6b0b2ab7f9e0b1ace309c72e767b2af071c4405e",