Last active
March 15, 2017 09:20
-
-
Save Cadair/5893586 to your computer and use it in GitHub Desktop.
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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# -*- 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()) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
An implementation of this is being discussed in sunpy/sunpy#1530