aiaprep returns wrong array size

[larrymanley]

I noticed that 'aiaprep' returns a 4098x4098 array rather than a 4096x4096 array.

[Cadair]

intresting, does map.rotate() do this?

[larrymanley]

I believe so; there's some padding in there, not sure if it's necessary...

[Cadair]

intersting. I will try and take a look. I guess this is a part of @ayshih change where he added the padding so the image is not cropped.

[wafels]

Perhaps we need a "keepdims" flag that returns an image the same size as
the original non-rotated image? For example, many LASCO science data
images taken in the SOHO rolled pointing era are rotated 7 degrees relative
to solar north. The LASCO processing routines rotate the data so that
solar north is "up" but the images are still 1024 x 1024, the same size
image as before SOHO had to roll. I suppose someone made the decision that
1k by 1k images were more important than the small amount of data that ends
up missing in the resulting rotated images.

On Tue, Mar 1, 2016 at 4:02 PM, Stuart Mumford notifications@github.com
wrote:

intersting. I will try and take a look. I guess this is a part of @ayshih
https://github.com/ayshih change where he added the padding so the
image is not cropped.

—
Reply to this email directly or view it on GitHub
#1691 (comment).

[ayshih]

I guess this is a part of @ayshih change where he added the padding so the image is not cropped.

I presume that is precisely the case. I'm not sure I like @wafels's suggestion of a keepdims flag, although I'm not strongly opposed to it either. It just seems to me that we shouldn't have rotate() do much more than rotating and preserving everything. If a processing routine such as aiaprep wants to enforce the same dimensions, it should explicitly submap for those dimensions as a separate statement following the rotation correction.

[larrymanley]

Is it one pixel on all four sides that is padded? As far as submapping, does that have implications for memory fragmentation when processing large amounts of data? AIA images have nothing but noise on the outer 1-2 pixels. How about a non-padded form of rotate internal to aiaprep? I would also recommend being able to specify 'missing' and 'order' directly from aiaprep.

[ayshih]

Is it one pixel on all four sides that is padded?

To be clear, it's not padding just for the sake of padding. The dimensions of the data array are expanded to ensure that no data is cropped out during the rotation. For example, for the "worst case" of a 45-degree rotation, the dimension would be expanded up by a multiplicative factor of sqrt(2) (so, here, to ~5793 pixels on a side). Because no data is cropped out, the expansion is effectively symmetric around the center of the map (which may not be the same as the origin of the coordinate system).

As far as submapping, does that have implications for memory fragmentation when processing large amounts of data?

There is no memory fragmentation considerations for submapping; you're creating a whole new Map object. There will be a small performance hit from the additional step of object creation, but that's unlikely to be a significant component compared to other processing (e.g., the rotation itself).

How about a non-padded form of rotate internal to aiaprep?

Definitely not a fan of this. I think that my suggestion is the most Pythonic approach.

I would also recommend being able to specify 'missing' and 'order' directly from aiaprep.

You should probably create a new issue so that your request doesn't get buried.

[Cadair]

[...]Because no data is cropped out, the expansion is effectively symmetric around the center of the map (which may not be the same as the origin of the coordinate system).

I thought map.rotate always rotated around the centre of the HPC system (i.e. centre of disk?)

How about a non-padded form of rotate internal to aiaprep?

Definitely not a fan of this. I think that my suggestion is the most Pythonic approach.

I agree with @ayshih, cropping in aiaprep is the easiest solution for maintanability.

I do wonder what is the most physically correct way to crop, is it to crop to 4096x4096 or is it to crop to some physical arcsecond size? I suppose if we are cropping after the scaling to 0.6" has taken place they are the same thing, we would be cropping to 4096*0.6"=2457.6" across?

[Cadair]

ping @drewleonard42

[SolarDrew]

I thought map.rotate always rotated around the centre of the HPC system (i.e. centre of disk?)

If I remember rightly, it does by default, but doesn't have to. I could be wrong, I've not looked at it in a while.

the expansion is effectively symmetric around the center of the map (which may not be the same as the origin of the coordinate system).

Presumably this means when we come to crop the image back again, we'll have to do so around the centre of the array rather than the origin?

[Cadair]

Presumably this means when we come to crop the image back again, we'll have to do so around the centre of the array rather than the origin?

Well aiaprep aligns the centre of the disk to the centre of the image no?

[SolarDrew]

Excellent point. I should not try to have these conversations before I've woken up properly....

[ayshih]

I thought map.rotate always rotated around the centre of the HPC system (i.e. centre of disk?)

What I meant was that the default behavior of Map.rotate() preserves all of the data, so the data array itself is necessarily is rotated about the center of the array. The definition of the coordinates after the rotation uses the reference coordinate in the metadata (which may or may not be the HPC origin) as the fixed point of the rotation.

Of course, as pointed out above, aiaprep() uses the non-default setting of recenter=True to shift the reference coordinate to the center of the data array, so the padding is actually no longer symmetric and some data will get clipped unavoidably.

I do wonder what is the most physically correct way to crop, is it to crop to 4096x4096 or is it to crop to some physical arcsecond size? I suppose if we are cropping after the scaling to 0.6" has taken place they are the same thing, we would be cropping to 4096*0.6"=2457.6" across?

In a related vein, my snarky response to this issue as written would have been "And...?". What do users actually need here? Do they need pixel dimensions of exactly 4096 on a side? Do they need data dimensions of a specific spatial extent?

[Cadair]

@ayshih as far as I see it the argument here goes that the rotation will never crop useful data off if the array size is maintained (there is no data in the corners where the padding saves the data), at least for full disk AIA images. So why bother padding?

[ayshih]

Agreed, for aiaprep(), there's little justification for preserving any data at the fringes, so certainly there can be an instrument-specific decision to crop after the rotation as needed.

This issue says that aiaprep() is supposed to return data arrays of exactly 4096x4096. I just wanted clarification if that is the actual requirement, given that part of the point of the Map class is provide an abstraction of the data-storage format.

[Cadair]

@ayshih I agree that it should be an instrument specific decision. Map.rotate is fine the way it is as far as I am concerned.

[larrymanley]

[larrymanley]

Can someone come up with a systematic way to submap 4096x4096 coming out of aiaprep? The above values will vary slightly depending on time of year and point in orbit, so the padding could potentially go a little above 4108.

[ayshih]

As @Cadair as much suggested, you can crop to the "central" 4096x4096 region (since recenter=True will ensure that the center of the data array is the HPC origin), or you can specify the appropriate angular extent in physical coordinates (which should be exact since the pixel angular size is rescaled to 0.6"). The former approach would requires a hair more math in aiaprep(), but is arguably preferred if the requirement is for particular pixel dimensions as opposed to arcsecond dimensions. (Even if the end result is the same, the code should reflect however the requirement is worded for future maintainability.)

[Cadair]

Thinking about this, would the correct logic be the follwing pesudo code:

centre = np.floor(CRPIX[0])
left_x = centre-4096/2
right_x = centre+4096/2

out = in[left_y:right_y, left_x:right_x]

?

[larrymanley]

That will work.

[larrymanley]

submap([left, right] * u.pix, [left, right] * u.pix)

[larrymanley]

Will anything in the tests need to be changed/added to implement this?

[Cadair]

@larrymanley just a new test for aiaprep to make sure it doesn't regress at some point in the future.

[Cadair]

@larrymanley why is AIA_171_ROLL_IMAGE not up to the task of being a test image?

[larrymanley]

@Cadair Can you give me some guidance on how to use the roll image in a test? It's in the sample data, not the test data. I tried downloading it directly as an online test, but it's a fits.gz file and wouldn't open. It worked a few weeks ago, and I know astropy will open .gz files.

[Cadair]

I would do something like the following:

@pytest.mark.online
def test_aiaprep():
    sunpy.data.download_sample_data(overwrite=False)
    aiamap = sunpy.map.Map(AIA_171_ROLL_IMAGE)