Instantly share code, notes, and snippets.

# Cadair/gist:5893586 Last active Mar 15, 2017

How to compensate for differential rotation in the default sunpy image. This is based upon a combination of pr/452 and pr/482 my branch for which can be found here: https://github.com/Cadair/sunpy/tree/map_refactor_diff_rot
 # -*- coding: utf-8 -*- """ Created on Sat Jun 29 21:31:58 2013 @author: stuart """ import numpy as np from datetime import timedelta import matplotlib.pyplot as plt plt.ion() from skimage import transform from skimage.util import img_as_float import sunpy import sunpy.wcs as wcs import sunpy.coords.util as cu from sunpy.time import parse_time data = sunpy.Map(sunpy.AIA_171_IMAGE) #Scikit-image needs float arrays between 0 and 1 def to_norm(arr): arr = img_as_float(arr) if arr.min() < 0: arr += arr.min() arr /= arr.max() return arr def un_norm(arr, ori): arr *= ori.max() if ori.min() < 0: arr -= ori.min() return arr #Define a function that returns a new list of coordinates for each input coord. def warp_sun2(xy, data): x,y = xy.T #Define some tuples scale = [data.scale['x'], data.scale['y']] ref_pix = [data.reference_pixel['x'], data.reference_pixel['y']] ref_coord = [data.reference_coordinate['x'], data.reference_coordinate['y']] #Calculate the hpc coords hpc_coords = wcs.convert_pixel_to_data(data.shape, scale, ref_pix, ref_coord) #Do the diff rot rotted = cu.rot_hpc(hpc_coords[1], hpc_coords[0], data.date, parse_time(data.date)- timedelta(days=2)) #Go back to pixel coords x2,y2 = wcs.convert_data_to_pixel(rotted[0], rotted[1], scale, ref_pix, ref_coord) #Restack the data to make it correct output form xy2 = np.column_stack([x2.flat,y2.flat]) #Remove NaNs mask = np.isnan(xy2) xy2[mask] = 0.0 return xy2 #Do the warp: out = transform.warp(to_norm(data.data), inverse_map=warp_sun2, map_args={'data':data}) #de norm out = un_norm(out, data.data) plt.figure() plt.subplot(121) plt.imshow(data.data, cmap=data.cmap, origin='lower', norm=data.norm()) plt.subplot(122) plt.imshow(out, cmap=data.cmap, origin='lower',norm=data.norm())

### dpshelio commented Aug 13, 2015

 An implementation of this is being discussed in sunpy/sunpy#1530