Create a gist now

Instantly share code, notes, and snippets.

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

What would you like to do?
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())

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment