astropy.wcs: access to NAXIS1 and NAXIS2

Hi,

I am using astropy.wcs and I think that NAXIS1, NAXIS2 are not easily accessible although they are used by WCS.calc_footprint(). See http://docs.astropy.org/en/stable/wcs/

Yes, I can access them via

NAXIS1= hdulist[ext].header['NAXIS1']

but it would seem generally useful to be able to also access them via the wcs package.

wcs.print_contents() does not seem to list NAXISn or axes although calc_footprint must be getting at them.

It would be convenient to be able to use something like:

    print('NAXIS: ', WCS.wcs.naxis) WORKS

    print('NAXES: ', WCS.wcs.naxes)

    print('NAXIS1: ', WCS.wcs.naxis1)
    print('NAXIS2: ', WCS.wcs.naxis2)

    print('CRPIX: ', WCS.wcs.crpix) WORKS
    print('CDELT: ', WCS.wcs.cdelt) WORKS

Just a quick note: NAXIS1 and NAXIS2 were made private in #323

@nden

cc @nden

Just a sidenote, this has come up previously in other packages, e.g.: https://github.com/aplpy/aplpy/pull/284/files
I was sure it had another issue on astropy, but I can't find one, so I guess this is new! Anyway, I agree that a WCS-like object that contains image size info would be useful.

And to expand the discussion a little further, a possible option is to provide a WCS subclass that knows about the image shape, e.g.:

class WCSFrame(WCS):
    def __init__(self,...):
        super().__init__()
        self.naxis1 = something
        self.naxis2 = something
        self.naxis3 = something
        # or just an abstract shape
        self.shape

which would be a light wrapper that would parse these keywords from a FITS header on instantiation and would provide the right sizes when sliced. We have a few use cases already: at least, the astropy one, but probably many others. For example, comparing WCSes, you may really care about the shape of the image also being the same

@astrofrog

As @astrofrog mentioned astropy.wcs changed a couple of times from exposing to hiding the shape of the image ( #323) . There are good reasons to go either way. A slightly different approach is taken in GWCS, where a domain is a list of dictionaries, where each item corresponds to an axis:

[{'lower': float, 'upper': float, 'includes_lower': bool, 'includes_upper': bool, 'step': float}]

Perhaps something like this can be useful in astropy.wcs

In my opinion NAXIS1,2 should be removed completely from the WCS including calc_footprint which should not belong to WCS. I think calc_footprint is probably the only place in WCS where NAXIS is needed. While WCS is closely tied with an underlying image [and I would have preferred this approach], it appears that PyWCS and now astropy.wcs consider WCS to be an independent object that can live all by itself in some abstract WCS space except for need for NAXIS. I think now it is the time to cut the umbilical cord and let WCS live by itself.

If calc_footprint() is useful, it could be placed in wcsutil and have a form such as calc_footprint(wcs, image_data).

A domain of validity of transformations in a WCS is a nice thing but it is not the same as the dimension/shape of the image's data array which is a property of the image and not WCS. For this reason I think it would be better to get rid of NAXIS1,2. Users should do ny, nx = image_data.shape and not expect WCS to do this for them as in the current astropy.wcs paradigm image data and WCS are independent things.

👎 to removing calc_footprint because I use it.

@pllim

@pllim stop using it! It's unhealthy :)

I need it to calculate footprint for visualization. Unless you have the equivalent functionality to replace this, I am opposed to removing it.

@pllim

@pllim Your visualization example is a perfect example when that shows that foot prints should be calculated from image AND WCS and not from WCS alone.

What about my proposal to move calc_footprint to wcsutils and have a function like this: calc_image_footprint(wcs, image_data)?

@mcara

@mcara I very much want access to an object that knows about the image shape and its WCS. It is important that the shape be an attribute of the same object, and it is important to me that I do not have to have the image (or a representation of the image) in memory. But I'm happy enough if the object that has this information is not the astropy.wcs.WCS object we have now, but some sibling class. Any thoughts on taking that approach?

Again, this is a personal opinion, but I think that fundamentally WCS is tied to an image (and its content - not only its shape). Therefore, your approach is not much different from what is it now. In fact, if we still want to mix WCS with an image, instead of introducing a new class I would rather see ._naxis1 being made public and have it printed to the header in the existing WCS.

Nowhere in the WCS standard (first one or two papers on WCS) could I find anything that would indicate that a WCS needs to know naxis1,2 to be properly defined. Thus, it is strange that anyone needs this info from a WCS and not from the image associated with that WCS.

@keflavich

@keflavich However, I can see your approach perfectly valid if you do it in your own code. That is, if my code was relying on WCS having naxis1,2 then I would have probably subclassed astropy.wcs.WCS just as you did in your example.

I disagree with the WCS being associated with an image's content. My most common use case for WCS is reprojecting one image into the exact pixel space of another, which means that the WCS is independent from the content.

But I also want to be very clear that I'm not making a semantic argument - I don't care if WCS is the object that knows about naxis + other WCS things, as long as there is some object that has both. I don't want to re-implement this multiple times, as would be required if I always subclass.

In the WCSFrame approach I suggested, one of the most important methods is to_header, which WCS also has, but WCSFrame.to_header() would create a fully-valid FITS header that specifies an image's shape too. This is what I need to programatically create a FITS header to reproject multiple data sets into. But this sort of object could also be passed directly to reproject instead of a fits.Header.

keflavich wrote: "My most common use case for WCS is reprojecting one image..." That is my point: you are re-projecting an image (not a WCS) and WCS is a tool to achieve this re-projection. Currently WCS provides you with everything you need to do this re-projection of an image even without NAXIS1,2 stored in the WCS. Size of the underlying image is not necessary to do everything you need otherwise it would have been part of the WCS standard.

It might be useful if I specified ONE of my use cases.

Often I want to determine the location of the boundaries of images to determine the image that an object lies on. The WCS header information for images or image extensions is stored in a metadata table in a database or in a fits catalogue.

It is therefore useful to be able to create a WCS header object using the header keywords. NAXIS1 and NAXIS2 are essential for this operation.

This similar to "Building a WCS structure programmatically": http://docs.astropy.org/en/stable/wcs/

I am coming for IDL where I could do this with:

http://idlastro.gsfc.nasa.gov/ftp/pro/astrom/make_astr.pro

@mcara

@mcara - what @keflavich meant is that the target WCS needs to have a shape/size - that is, currently in reproject if you specify a target WCS, you also need to specify a target shape, so he just wants a convenience class that contains WCS + shape (equivalent to a FITS header)

@astrofrog

@astrofrog @keflavich: "...that is, currently in reproject if you specify a target WCS, you also need to specify a target shape..." I did not know what @keflavich is working on so I may have misunderstood some parts.

However, unless we really tie together an image with its WCS, I would prefer to have them completely "independent" (at least formally). Specifically, let astropy.wcs be a "World Coordinate System" - that is, a set of coordinate transformations as it was intended to be. The image size to which this transformation is going to be applied has nothing to do with the coordinate transformations. Whether the image is 1x1 or 100,000x100,000 - it's all the same (from the coordinate transformation point of view). At the very least, conceptually, image size has nothing to do with coordinate transformation.

I understand that having naxis1,2 stored with the WCS may be convenient in some cases but I do not think it is necessary. A long time ago I also was relying on WCS having naxis1,2 in my code. So I was pretty upset when naxis was made private. However, I have moved on and I do not even care what naxis1,2 is in a WCS object.

The problem with having image properties mixed with coordinate transformations is that things become a mess and lead to issues such as #4498 or someone could cut the image, or rescale it, etc. and not update naxis1,2 in the WCS. It is my belief that things are much cleaner when we conceptually separate image data/properties from coordinate transformations (unless we tightly unify these concepts in a single class ImageOnTheSky).

@keflavich

Regarding @keflavich idea of subclassing WCS (class WCSFrame(WCS))... I think it would be fine to do this in your project and require users to explicitly "tie" a WCS with image properties (shape). I just would not use "Frame" in the name as "frame" hints (to me) to a "coordinate frame" (detector, sky, etc.) I have done something similar to deal with limitations of old pyregion not using distortions defined in WCS - see class _AuxSTWCS here: https://github.com/spacetelescope/drizzlepac/blob/master/lib/drizzlepac/mapreg.py

@richardgmcmahon

@richardgmcmahon Let me address some issues mentioned by you in a previous post. You said:

It is therefore useful to be able to create a WCS header object using the header keywords. NAXIS1 and NAXIS2 are essential for this operation.

This similar to "Building a WCS structure programmatically": http://docs.astropy.org/en/stable/wcs/

Part 1.

NAXIS1 and NAXIS2 is required to create a FITS header containing data and so NAXIS1,2 should indicate the shape of data. If you want to save a WCS to a fits file without data - you can still do it as you will see soon in the next part.

Part 2.

Let's go through the example in the "Building a WCS structure programmatically" link (and let's also play with that WCS's footprint - just for fun).

>>> import numpy
>>> from astropy import wcs
>>> from astropy.io import fits
>>> w = wcs.WCS(naxis=2)
>>> w.wcs.crpix = [-234.75, 8.3393]
>>> w.wcs.cdelt = numpy.array([-0.066667, 0.066667])
>>> w.wcs.crval = [0, -90]
>>> w.wcs.ctype = ["RA---AIR", "DEC--AIR"]
>>> w.wcs.set_pv([(2, 1, 45.0)])
>>> pixcrd = numpy.array([[0, 0], [24, 38], [45, 98]], numpy.float_)
>>> world = w.wcs_pix2world(pixcrd, 1)
>>> print(world)
[[ 267.96547027  -73.73660749]
 [ 276.53931377  -71.97412809]
 [ 287.77080792  -69.67813884]]
>>> pixcrd2 = w.wcs_world2pix(world, 1)
>>> print(pixcrd2)
[[  1.16529009e-12  -2.80664381e-13]
 [  2.40000000e+01   3.80000000e+01]
 [  4.50000000e+01   9.80000000e+01]]
>>> header = w.to_header()
>>> header
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =              -234.75 / Pixel coordinate of reference point            
CRPIX2  =               8.3393 / Pixel coordinate of reference point            
CDELT1  =            -0.066667 / [deg] Coordinate increment at reference point  
CDELT2  =             0.066667 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---AIR'           / Right ascension, Airy's zenithal projection    
CTYPE2  = 'DEC--AIR'           / Declination, Airy's zenithal projection        
CRVAL1  =                  0.0 / [deg] Coordinate value at reference point      
CRVAL2  =                -90.0 / [deg] Coordinate value at reference point      
PV2_1   =                 45.0 / AIR projection parameter                       
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =                -90.0 / [deg] Native latitude of celestial pole        
RADESYS = 'ICRS'               / Equatorial coordinate system                   
>>> hdu = fits.PrimaryHDU(header=header)
>>> hdu.header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                    8 / array data type                                
NAXIS   =                    0 / number of array dimensions                     
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =              -234.75 / Pixel coordinate of reference point            
CRPIX2  =               8.3393 / Pixel coordinate of reference point            
CDELT1  =            -0.066667 / [deg] Coordinate increment at reference point  
CDELT2  =             0.066667 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---AIR'           / Right ascension, Airy's zenithal projection    
CTYPE2  = 'DEC--AIR'           / Declination, Airy's zenithal projection        
CRVAL1  =                  0.0 / [deg] Coordinate value at reference point      
CRVAL2  =                -90.0 / [deg] Coordinate value at reference point      
PV2_1   =                 45.0 / AIR projection parameter                       
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =                -90.0 / [deg] Native latitude of celestial pole        
RADESYS = 'ICRS'               / Equatorial coordinate system

Please notice that:

NAXIS   =                    0 / number of array dimensions                     
WCSAXES =                    2 / Number of coordinate axes

that is, the image has 0 axes while the WCS has two axes. WCS does not need to know how "large" is each axis. In addition, notice that the header created with fits.PrimaryHDU(header=header) does not contain NAXIS1 or NAXIS2 !!! In order to get "reasonable" NAXIS1,2 values, you will need to provide "data" to fits.PrimaryHDU.

Finally, let's see what calc_footprints gives us with the above WCS:

>>> w.calc_footprint()
array([[ 268.21686048,  -73.66991904],
       [ 267.97409309,  -73.66762507],
       [ 267.96547027,  -73.73660749],
       [ 268.2092695 ,  -73.73891131]])
>>> w.all_world2pix(w.calc_footprint(),1)
array([[  1.00000000e+00,   1.00000000e+00],
       [  1.00000000e+00,  -1.54543045e-13],
       [  1.16529009e-12,  -2.80664381e-13],
       [  5.11590770e-13,   1.00000000e+00]])

So, the footprint of this image is 1x1 pixel large (should have been 0, in my opinion) - probably not even close to the size of your real image.

Just to add to my previous post. In your IDL script you are using information about BOTH coordinate transformation (WCS) AND image properties (in particular, image size) to produce a header that has NAXIS1,2. However, this header may be wrong if you do not save image data to the FITS file. The way it is currently done in astropy you do not need to worry about this if you do the things using, e.g., fits.PrimaryHDU():

>>> d = numpy.zeros((15,15))
>>> hdu = fits.PrimaryHDU(header=header, data=d)
>>> hdu.header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -64 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                   15                                                  
NAXIS2  =                   15                                                  
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =              -234.75 / Pixel coordinate of reference point            
CRPIX2  =               8.3393 / Pixel coordinate of reference point            
CDELT1  =            -0.066667 / [deg] Coordinate increment at reference point  
CDELT2  =             0.066667 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---AIR'           / Right ascension, Airy's zenithal projection    
CTYPE2  = 'DEC--AIR'           / Declination, Airy's zenithal projection        
CRVAL1  =                  0.0 / [deg] Coordinate value at reference point      
CRVAL2  =                -90.0 / [deg] Coordinate value at reference point      
PV2_1   =                 45.0 / AIR projection parameter                       
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =                -90.0 / [deg] Native latitude of celestial pole        
RADESYS = 'ICRS'               / Equatorial coordinate system

Notice that PrimaryHDU header now has NAXIS=2 and NAXIS1,2=15. So, it is PrimaryHDU that takes care of putting correct NAXIS1,2 based on image data.

Thanks for the clarification;

just to be clear; for this use case; i do not write an image.

I am building a WCS from headers keywords that are within a table and using WCS to determine the footprint of a image. I take each RA, Dec that I am working with;
(i) do a radial pairwise search on the image metadata using an inscribed radius defined by the image corners.
(ii) then use WCS to compute the pixel a position has on an image; test against the image shape as defined by NAXIS1, NAXIS2

It would convenient to use the use the WCS object as a container that includes he NAXIS1, 2 since I can then pass the WCS and
NAXIS1, 2 info as single argument.

I understand the rationale for separation of image shape from WCS transformation and do not have anything significant to contribute.

I welcome very much the very comprensive support for WCS in astropy as someone who deals with imaging data from a multitude of imaging surveys and wanted to move away from IDL and bespoke code.

I also understand that the concept of image WCS and image properties are different concepts and also the risk of creating images with the
wrong headers when one trims or cutouts out a subset.

Hopefully my use case has some merit in the planning of features. I am not blocked at the moment.

regards Richard

On 20 Apr 2016, at 00:01, Mihai Cara [email protected] wrote:

Just to add to my previous post. In your IDL script you are using information about BOTH coordinate transformation (WCS) AND image properties (in particular, image size) to produce a header that has NAXIS1,2. However, this header may be wrong if you do not save image data to the FITS file. The way it is currently done in astropy you do not need to worry about this if you do the things using, e.g., fits.PrimaryHDU():

d = numpy.zeros((15,15
))

hdu = fits.PrimaryHDU(header=header, data=
d)

hdu.header

SIMPLE = T / conforms to FITS
standard

BITPIX = -64 / array data type

NAXIS = 2 /
number of array dimensions

NAXIS1 = 15

NAXIS2 = 15

WCSAXES = 2 /
Number of coordinate axes

CRPIX1 = -234.75 /
Pixel coordinate of reference point

CRPIX2 = 8.3393 /
Pixel coordinate of reference point

CDELT1 = -0.066667 /
[deg] Coordinate increment at reference point

CDELT2 = 0.066667 /
[deg] Coordinate increment at reference point

CUNIT1 = 'deg' / Units of coordinate increment and
value

CUNIT2 = 'deg' / Units of coordinate increment and
value

CTYPE1 = 'RA---AIR' / Right ascension, Airy's zenithal projection
CTYPE2 = 'DEC--AIR' / Declination, Airy's zenithal projection
CRVAL1 = 0.0 /
[deg] Coordinate value at reference point

CRVAL2 = -90.0 /
[deg] Coordinate value at reference point

PV2_1 = 45.0 / AIR
projection parameter

LONPOLE = 180.0 /
[deg] Native longitude of celestial pole

LATPOLE = -90.0 /
[deg] Native latitude of celestial pole

RADESYS = 'ICRS' / Equatorial coordinate system
Notice that PrimaryHDU header now has NAXIS=2 and NAXIS1,2=15. So, it is PrimaryHDU that takes care of putting correct NAXIS1,2 based on image data.

—
You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub

@mcara

@mcara - I think that there are genuine use cases where one wants to represent WCS + an image extent without wanting to set specific pixel data or create a FITS header (which is FITS-specific). So to move this discussion forward, I'd like to suggest we start thinking about @keflavich's suggestion to make a WCS subclass. I understand (and agree that) WCSFrame may not be a good name, so does anyone have any other ideas of what this could be? We could then remove entirely any concept of size from the WCS class (including in private attributes).

I believe this has actually been addressed as pixel_shape and array_shape are exposed and optional in APE 14:

In [1]: from astropy.wcs import WCS

In [2]: wcs = WCS('filename.fits')

In [3]: wcs.array_shape
Out[3]: (2000, 1000)

In [4]: wcs.pixel_shape
Out[4]: (1000, 2000)

Navigation Menu

Repository navigation

astropy.wcs: access to NAXIS1 and NAXIS2 #4669

Part 1.

Part 2.

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

astropy.wcs: access to NAXIS1 and NAXIS2 #4669

Description

Activity

astrofrog commented on Mar 5, 2016

astrofrog commented on Mar 5, 2016

keflavich commented on Mar 27, 2016

keflavich commented on Mar 28, 2016

nden commented on Apr 19, 2016

mcara commented on Apr 19, 2016

mcara commented on Apr 19, 2016

pllim commented on Apr 19, 2016

mcara commented on Apr 19, 2016

pllim commented on Apr 19, 2016

mcara commented on Apr 19, 2016

keflavich commented on Apr 19, 2016

mcara commented on Apr 19, 2016

mcara commented on Apr 19, 2016

keflavich commented on Apr 19, 2016

mcara commented on Apr 19, 2016

richardgmcmahon commented on Apr 19, 2016

astrofrog commented on Apr 19, 2016

mcara commented on Apr 19, 2016

mcara commented on Apr 19, 2016

mcara commented on Apr 19, 2016

Part 1.

Part 2.

mcara commented on Apr 19, 2016

richardgmcmahon commented on Apr 20, 2016

astrofrog commented on Apr 20, 2016

astrofrog commented on Sep 23, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions