Independent thresholding for each layer in composite maps

[alexjsolar]

When using the composite map feature, I'd like to be able to set different vmin and vmax constraints for each layer. At the moment, I can set an overall vmin and vmax for the plot, but the scales of the plots used in each layer are different, so one-size-fits-all thresholds aren't appropriate.

Something similar to the alpha value that can be set for each layer would be great.

[dpshelio]

I think that adding somthing like set_minmax(layer, vmin, vmax) or similar could work.
set_levels() does something similar.

[abit2]

Hey i would like to give it a try

[abit2]

@ajamesucl @dpshelio so need to set different contour threshold levels in each layer and if it is less than vmin than =vmin else equal to level. Similar for vmax .
Did i understand it right please correct me if wrong

[alexjsolar]

Just to be clearer about what I specifically would like to do: in my example, the composite contains two layers. The first being a map of the Sun, and then contours taken from another map are overlayed as a second layer. For example, to plot magnetic field strength contours from an HMI magnetogram on top of an AIA image of the Sun.

The problem is that I would like to plot each layer with different thresholds. So I would like to plot between a 'vmin1' and 'vmax1' for the base layer, and plot the second layer between a 'vmin2' and a 'vmax2'.

But yes I think you understand right about what the thresholding does. Thanks!

[abit2]

def set_minmax(self, index, levels, vmin, vmax):
"""
Sets the contour threshold levels for a layer in the composite image.

    Parameters
    ----------
    index : `int`
        The index of the map in the composite map.

    levels : array-like
        The contour levels.

    vmin :   
        Minimum contour threshold level.

    vmax :
        Maximum contour threshold level.

    Returns
    -------
    `~sunpy.map.CompositeMap`
        A composite map with contour levels 'levels' at layer 'index'.
    """
    for level in levels:
        if level > vmax:
            self._maps[index].levels.append(vmax)
        if level < vmin:
            self._maps[index].levels.append(vmin)
        else:
            self._maps[index].levels.append(level) 

is this right?? @ajamesucl
Please correct if wrong

[dpshelio]

@abit2 you should do this as a Pull Request and link in the comment with this issue number, that way the repository is tested automatically (and checks whether that breaks something) and we can provide inline comments.

I think you are confusing the two things composite maps can do, one is to overlay contours on top of an image (controlled by levels) or you can visualise two images with semi-transparency (alpha value). The case @ajamesucl is talking about is the last one, and does not have to do with levels. The way I would do it is similar as set for alpha - as that is a matplotlib parameter.

[larrymanley]

@ajamesucl Perhaps you are aware of this, but if you are mixing AIA and HMI data, the two instruments have different plate scales (CDELT*) - AIA is ~0.6 arcsec/pixel, HMI is ~0.5, so the HMI data needs to be scaled down. I'll be doing an "hmiprep" with a plate scale option and optional removal of the limb noise on magnetograms.

[alexjsolar]

@larrymanley Thanks for the help! I just checked the CDELTs for my AIA and HMI images and both are the same at ~0.6 arcsec/pixel. I've been getting and prepping my data with IDL Solarsoft, so I had used the aiaprep function on both sets of data, which must have accounted for the HMI scale.

[sudk1896]

@ajamesucl : I want to take a stab at this. Could you give me a use case, specifically, any two fits files for the maps that you're using for this feature (#1697 (comment)). This would help me understand and work out the plot more clearly.
Thanks.

[alexjsolar]

Hi @sudk1896 I've included a Dropbox link to an HMI fits file and an AIA fits file that I'd like to turn into a composite map, as well as some basic code to do this: https://www.dropbox.com/sh/ws1e2bsr9nu0p9s/AADUIu0dDy7bLjc6OG2s6XsTa?dl=0

My request is that I would like to plot each map in the composite between different limits. For example, I would like to plot the HMI map layer between vmin1=-2000 and vmax1=+2000, and plot the AIA map layer on top of this between vmin2=0 and vmax2=1500.
Thanks for your help!

[sudk1896]

@ajamesucl: No problem. I would get started on this.

[larrymanley]

Is this close to what you have in mind? I need to make a few refinements and I'll send you the code.


This image shows the difference in plate scale (I used your Dropbox files):

I'm not sure why some of the AIA data is plotting in grey.

[alexjsolar]

That does look like what I was hoping for! Thanks!

I tried changing the scale of that HMI image from 0.5 arcsec/pixel to 0.6 so it matched AIA, but this actually seemed to throw the alignment of features off.

[sudk1896]

@ajamesucl: Do you still want me to work on this ?

[dpshelio]

The code needs to be implemented, @larrymanley have you changed how map behaves? would you be able to do a Pull Request with the changes?
if not, @sudk1896 could make it for a general case.

[larrymanley]

A few notes:
I used an inline equivalent of aiaprep for the HMI image:

hmi_scale_factor = hmimap.scale.x / (0.6 * u.arcsec)
hmimap = hmimap.rotate(recenter=True, scale=hmi_scale_factor.value, missing=0.0)

The first two images above are still off a little bit because of the padding in rotate (issue #1691). I've got a fix for aiaprep but I need a full-res (4096x4096) level 1 file for the the test.

I used Numpy clip on both maps rather than using vmin/vmax in the plot settings:

aiamap.data = np.clip(aiamap.data, 0, 1500)
hmimap.data = np.clip(hmimap.data, -2000.0, 2000.0)

When I do my AIA movies I pop the 'norm' out of plot_settings because I do "LMSAL-style" intensity scaling. It doesn't appear that compositemap works without the 'norm' in the plot_settings of the AIA image. I think this type of composite is a good example of why 'norm' should be optional (line 460 of compositemap.py).

[dpshelio]

Thanks @larrymanley for sharing!! Thanks for the norm comment, I agree that it should be optional. Maybe we should add another issue to don't forget about this.

@sudk1896 what I would like to see implemented is within the map object, so the user doesn't have to clip the data before plotting, but change the clipping at plotting time. It's not very difficult, it's similar to set_levels().

[Cadair]

norm should just default to None, @larrymanley if you want no norm, then set it to None I think that will do what you expect?

[larrymanley]

@ajamesucl Here's an example script. I'll experiment with vmin/vmax plot_settings tonight.

import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
import sunpy.map
from sunpy.instr.aia import aiaprep

hmimap = sunpy.map.Map('hmi.fits')
aiamap = sunpy.map.Map('aia.fits')

hmi_scale_factor = hmimap.scale.x / (0.6 * u.arcsec)

hmimap = hmimap.rotate(recenter=True, scale=hmi_scale_factor.value, missing=0.0)
# extract center 4096x4096 array
center = np.floor(hmimap.meta['crpix1'])
leftx = center - 4096/2
rightx = center + 4096/2
hmimap = hmimap.submap([leftx, rightx]*u.pix, [leftx, rightx]*u.pix)
hmimap.meta['r_sun'] = hmimap.meta['rsun_obs'] / hmimap.meta['cdelt1']

hmimap.data = np.clip(hmimap.data, -2000.0, 2000.0)

# remove HMI magnetogram off-limb noise
# thanks Drew Leonard!
lowx, highx = (hmimap.xrange[0].value / hmimap.scale.x.value,
               hmimap.xrange[1].value / hmimap.scale.x.value)
lowy, highy = (hmimap.yrange[0].value / hmimap.scale.y.value,
               hmimap.yrange[1].value / hmimap.scale.y.value)
x_grid, y_grid = np.mgrid[lowx:highx-1, lowy:highy-1]
r_grid = np.sqrt((x_grid ** 2.0) + (y_grid ** 2.0))
outer_rad = (hmimap.rsun_obs.value * 1.0) / hmimap.scale.x.value
hmimap.data[r_grid > outer_rad] = hmimap.min()

hmimap.peek(draw_limb=True)

aiamap = aiaprep(aiamap)
# aiamap.data = np.clip(aiamap.data, 0.0, 1500.0)

aiamap.data = np.sqrt(np.clip((aiamap.data*(4.99803 / aiamap.exposure_time.value)), 10, 6000))
# LMSAL/The Sun Today intensity scaling
aiamap.plot_settings['norm'] = None

cmap1 = hmimap.plot_settings['cmap'] = plt.get_cmap('Greys_r')
cmap2 = aiamap.plot_settings['cmap'] = plt.get_cmap('sdoaia171')

comp_map = sunpy.map.Map(hmimap, composite=True, alpha=1.0)
comp_map.add_map(aiamap, alpha=0.25)

fig = plt.figure()
ax = fig.add_subplot(111)
comp_map.plot()
plt.colorbar()
plt.show()

[sudk1896]

@ajamesucl: Can you post the dropbox link again. It has expired. I would like to resume where I left off. Thanks.

[dpshelio]

@sudk1896 you don't need to use his data for this, you can use the sample data in SunPy, take any two maps, if you want to get closer to his example use the AIA171 and the HMI

[Cadair]

@wafels is this closed by #1906 ?

[wafels]

Yes, #1906 closes this. This is because there are now get and set methods
for the plot_settings attribute for each map in the composite, and vmin and
vmax can be set through the image normalization method held in
plot_settings.

On Thu, Oct 13, 2016 at 7:33 AM, Stuart Mumford notifications@github.com
wrote:

@wafels https://github.com/wafels is this closed by #1906
#1906 ?

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
#1697 (comment), or mute
the thread
https://github.com/notifications/unsubscribe-auth/AA8CF7KFM8A3vFwWMpigsCC2IOT9MxXFks5qzhcMgaJpZM4HoVjC
.