Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Lack of rotation when plotting rotated fields of view with NDCube #99

Closed
BaptistePellorceAstro opened this issue Mar 22, 2018 · 16 comments

Comments

@BaptistePellorceAstro
Copy link
Contributor

By testing some data to check if coordinates are well represented with NDCube instance, we found that for rotated fields of view, X and Y axis are not representing usual coordinates that we use. We tried with iris_l2_20150304_160141_3824263489_SJI_1400_t000 data coming from a 45 degree's rotated field of view, and you can see, below, the result we have. In this case, the images should be rotated by 45 degrees counter-clockwise.
iris_l2_20150304_160141_3824263489_SJI_1400_t000.fits

@BaptistePellorceAstro BaptistePellorceAstro changed the title Lack of rotation when plotting rotated fields of view Lack of rotation when plotting rotated fields of view with NDCube Mar 22, 2018
@Cadair
Copy link
Member

Cadair commented Mar 22, 2018

If you look closely the ticks are rotated, so the plot knows the coordinate system does not align to the image.

What needs changing here is the defaults for the tick labels and probably a faint grid to make it more obvious that the image and the coordinate system are not aligned.

The plot will not and should not rotate the image to the coordinate system.

@Cadair
Copy link
Member

Cadair commented Mar 22, 2018

Having said that, work is currently being done on astropy to make WCSAxes default to standard decimal notation for arcsec by default. (astropy/astropy#7318)

@tiagopereira
Copy link
Member

@Cadair so this is really an issue with WCSAxes and not ndcube? Maybe this is not a problem for astronomical images, but for solar observations I think it would be really useful to have the option of rotating the image so that the coordinates are aligned with the figure axes. It is very common to work with multi-telescope observations, and superimposing two observations on the same coordinate system is really useful.

@Cadair
Copy link
Member

Cadair commented Mar 22, 2018

@tiagopereira overlaying non-aligned data is not an operation that WCSAxes / matplotlib can or should just handle, it's one for a package like reproject which you can use to do that for display purposes with solar data.

I think there is an NDCube issue here, in that we should have some more sane defaults for the tick labels and grids etc, i.e. match the default settings of Map.plot in sunpy core.

@tiagopereira
Copy link
Member

@Cadair thanks for the link to reproject, which seems to be very neat (although I worry about all the package fragmentation: astropy / reproject / sunpy / irispy / ndcube). However, reproject re-grids the images, which is not exactly what I had in mind. In particular, for visualising with matplotlib some re-gridding/interpolation would of course be necessary, but ideally we would still preserve the original, non-interpolated data for any other operations. One could for example plot an AIA image and superimpose an IRIS slit-jaw image, keeping their different pixel scales and not necessarily scaling up or down to match pixel-by-pixel.

Of relevance to NDCube is that it appears that when the coordinates are not aligned with the image, its slicing crop_by_coords() method doesn't make a lot of sense. For example, using the IRIS file from Baptiste's plot:

a = fits.open("iris_l2_20150304_160141_3824263489_SJI_1400_t000.fits")
ww = WCS(a[0].header)
data = NDCube(a[0].data, ww)
amin = Quantity(1., unit='arcmin')
roi = data.crop_by_coords([Quantity(13063, unit='s'), amin*4, -amin*7], 
                          [Quantity(500, unit='s'), amin, amin])
roi.plot()

roi

@DanRyanIrish
Copy link
Member

What about crop_by_coords doesn't make sense to you, @tiagopereira? If it's the fact that the plot is not very viewable it's because plot() is using an aspect ratio equal to the number of pixels in each dimension. It's not to do with crop_by_coords. You could try doing
roi.plot(aspect="auto")
and that should expand the y-axis for you. If however, there's something else about crop_by_coords that I'm missing that doesn't make sense, could you elaborate to help me understand? Thanks!

Regarding the point about sane defaults for the tick labels and grids, ndcube currently has no specific defaults! We just take whatever matplotlib give us back as defaults. Improving this would definitely be a useful enhancement.

@tiagopereira
Copy link
Member

@DanRyanIrish the aspect ratio is just the symptom, the problem seems to be that crop_by_coords is trying to crop the minimum region in pixel space that satisfies the criteria, and in this extreme case (45 deg) it falls well short of the region of interest.

What I'm doing in the commands above is trying to crop a similar image to what Baptiste plotted. If you see his figure, the x axis goes from ~-7' to -5', and y from ~3' to 4.5'. In my roi, the requested crop was from -7' to -6' in x and 4' to 5' in y. This should give a diagonal slice off that image, with some white space where there is no data. Instead, what I get is a 6x510 image that covers the x range -7' to -6' (ok so far), but only 4' to ~4.01' on the top left corner. This "slice" was because the top right corner is at 5', thereby fulfilling the criteria for crop_by_coords. It would be nice to obtain the whole region satisfying those criteria, not just the minimum region. But ideally, what would be even better was to obtain some kind of boolean map that followed the coordinates on the original image, which eventually could be worked with reproject or subject to further processing.

@Cadair
Copy link
Member

Cadair commented Mar 22, 2018

@tiagopereira You do appear to have found a true bug in crop_by_coords I suggest we move that discussion to another issue.

As for the package fragmentation, that is a design decision for Python, it's a very modular ecosystem. Given it's very easy to install packages and have packages depend on one another, I really think it's a source of strength once you are used to it.

Overplotting images with different pixel scales I think is possible using matplotlib, but you need to re-grid them to align them before you plot them, because matplotlib does all the drawing in pixel space. We should look into wrapping this functionality up into IRISPy or sunpy where appropriate.

@DanRyanIrish
Copy link
Member

Ah, I better understand now. Thanks @tiagopereira. This reveals that crop_by_coords is best suited to cases where the pixel grid is aligned with the real world coordinate grid. It clearly could do with some more development for more complicated cases. I think your idea of getting the maximum region that satisfies the criteria (rather than the minimum) and using the mask to pull out non-rectangular regions of interest is a good one. I'm not sure crop_by_coords should do reprojection, but perhaps there's scope for another method that can do that, or perhaps a reproject=True kwarg? Although perhaps it's wiser to use reproject at this stage. I don't reprojection is a trivial issue.

Perhaps yourself or @BaptistePellorceAstro would be interested in contributing to ndcube to change crop_by_coords to your first suggested implementation as it's the simplest to start with, i.e. maximum non-projected region + masking.

Any thoughts @Cadair?

@Cadair
Copy link
Member

Cadair commented Mar 22, 2018

I don't think ndcube should use reprojection, that is a bigger picture problem that needs looking at in sunpy. As much as I do like the idea of having a reprojection aware overplotting helper in SunPy.

See issue #100 for discussion on the crop_by_coords issue. We should leave this one to the plotting labelling issues.

@DanRyanIrish
Copy link
Member

-- "I don't think ndcube should use reprojection, that is a bigger picture problem that needs looking at in sunpy. As much as I do like the idea of having a reprojection aware overplotting helper in SunPy."

Agreed!

@tiagopereira
Copy link
Member

Agree that allowing for general WCS coordinates, reprojection can be quite complex. But most of the cases people deal with in solar observations should not be that complicated -- it is just a rotation away. E.g., most space data comes from Earth orbiting satellites where the reprojection is simple. Many (most?) ground-based telescopes of the 1m+ class have alt-az mounts and no derotator, meaning that the solar image rotates with time in the CCD plane. If we also want to use NDCube for those data (I think we should), then this rotation is necessary for every frame.

This is going beyond the scope of the original discussion, but I think it is important to make a distinction between this rotation of the coordinate space for plotting and cropping purposes vs for storing and manipulating the data. Not only can the rotation operation (interpolation) be destructive on the data, storing rotated data can result in much larger file sizes because the array sizes increase (factor of 2 for 45 degrees, much more when you consider IRIS SJIs that cover a wide raster). For "plotting and cropping" I feel this is very much on the realm of NDCube.

@DanRyanIrish
Copy link
Member

I think that plotting and cropping are distinctly different things. Cropping is certainly within the scope of ndcube. As for plotting, I am currently rewritting NDCubeSequence so that the plotting is a mixin class, i.e. an optional set of extra functionalities. This is already the case for NDCube. This will mean that people --- solar and non-solar, space-based and ground-based --- will be able to write their own plotting mixins and combine them with ndcube depending on their needs. But we are currently planning to provide a ready-to-go mixin as the lines of what already exists.

I completely agree on with @tiagopereira on the usefulness of rotating for plotting purposes. I think it's a question of where that core functionality lives. I think due to its complexity, and because work has already been done on this topic for sunpy.map, such functionality should be developed in sunpy and then imported into the ndcube plotting mixins. This makes very little difference in the development of this functionality (it's just a matter of which repo you PR to), especially since the ndcube and sunpy community overlap so much. But it does make more sense for maintenance in the future.

My discussion here really has gone down the rabbit hole to the point of being off topic. My main point is that being able to access this functionality through NDCube (or a subclass of it) would be highly useful, as @tiagopereira pointed out. It's just that most the actual code that makes it possible should live in sunpy and then be imported in ndcube or a sublclass of one of NDCube or NDCubeSequence.

@Cadair
Copy link
Member

Cadair commented Mar 22, 2018

[I wrote this before Danny's comment]

So, in general, sunpy (core) has tried to work on the principal of doing everything in coordinate space without having to rotate (reproject) the data. And yes, while I realise that rotation and arbitrary reprojection differ in complexity, they are effectively the same underlying operation, so for the most part they are equivalent in this discussion.

[side note, for ndcube it would be a lot easier to just use the reproject package than trying to convert GenericMap.rotate to working on higher dimensional data]

I don't think you need reprojection for cropping in physical coordinates, you just need to make sure you take the maximal coordinate box as discussed here and in #100. (As for generating masks of coordinate regions, this is super easy using axis_world_coords.)

As for plotting, I think it would be quite easy to write a overplotting script for any NDData subclass (so NDCube and GenericMap) which uses WCSAxes and reproject. fwiw it would be worth looking at this WIP glue PR which is for using reproject to do automatic overplotting.

@DanRyanIrish
Copy link
Member

I'd just like to thank @BaptistePellorceAstro for raising this issue and @tiagopereira for his contributions to it. As we can see it's sparked some very interesting discussion about ndcube development as well as highlighting the bug outlined in Issue #100. This is just the sort of vibrant discussion and feedback ndcube can benefit from.

@tiagopereira
Copy link
Member

@Cadair sorry, I seem to have missed your earlier post about package fragmentation and plotting. I happen to be responsible for maintaining the python distributions in our group, and from that perspective it is often annoying to have many packages to install/update. Also I get similar feedback from new users that are trying to install in their own hardware. But anyway...

About matplotlib plotting in pixel scale, this is not always the case. For example imshow has the extent option, which allows arbitrary stretches in x and y (and effectively rectangular pixels), and pcolor/pcolormesh support varying grids. Perhaps more relevant to this example is mpl_toolkits.axisartist.floating_axes, which allows a simple rotation of the whole axes but at a loss of funcionality (e.g. zooming in and interactive visualisation). So it is possible in principle, but may have too many tradeoffs...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants