Dimensions and World Coordinates

AstroImages are based on Dimensional Data. Each axis is assigned a dimension name and the indices are tracked.

World Coordinates

FITS files with world coordinate system (WCS) headers contain all the information necessary to map a pixel location into celestial coordinates & back.

Let's see how this works with a 2D image with RA & DEC coordinates.

using AstroImages
using Plots

# Download a Hubble image of the Eagle nebula
download(
    "http://www.astro.uvic.ca/~wthompson/astroimages/fits/656nmos.fits",
    "eagle-656nmos.fits"
);
eagle = load("eagle-656nmos.fits")

This image contains world coordinate system headers. AstroImages.jl uses WCS.jl (and wcslib under the hood) to parse these headers. We can generate a WCSTransform object to inspect:

wcs(eagle, 1) # specify which coordinate system

WCSTransform(naxis=2,cdelt=[1.0, 1.0],crval=[274.71149247724, -13.816384007184],crpix=[386.5, 396.0])

Note that we specify with an index which coordinate system we'd like to use. Most images just contain one, but some contain multiple systems.

We can lookup a coordinate from the image:

world = pix_to_world(eagle, [1, 1]) # Bottom left corner

2-element Vector{Float64}:
 274.712299241082
 -13.801135972688115

Or convert back from world coordinates to pixel coordinates: We can lookup a coordinate from the image:

world_to_pix(eagle, world) # Bottom left corner

2-element Vector{Float64}:
 1.000000000336172
 0.9999999992196535

These pixel coordinates do not necessarily have to lie within the bounds of the original image, and in general lie at a fractional pixel position.

If an image contains WCS headers, we can visualize them using implot:

implot(eagle)

We can adjust the color of the grid:

implot(eagle, gridcolor=:cyan)

If these aren't desired, we can turn off the grid or the WCS tick marks:

plot(
  implot(eagle, grid=false),
  implot(eagle, wcsticks=false),
  size=(900,300),
  bottommargin=10Plots.mm
)

Since AstroImages are based on DimensionalData's AbstractDimArray, the mapping between pixel coordinates and world coordinates are preserved when slicing an AstroImage:

slice1 = eagle[1:800,1:800]
slice2 = eagle[800:1600,1:800]
plot(
  implot(slice1),
  implot(slice2),
  size=(900,300),
  bottommargin=10Plots.mm
)

World coordinate queries from that slice are aware of their position in the parent image:

@show pix_to_world(slice1, [1,1])

2-element Vector{Float64}:
 274.712299241082
 -13.801135972688115

@show pix_to_world(slice2, [1,1])

2-element Vector{Float64}:
 274.7277517847315
 -13.817350009028138

Note that you can query the dimensions of an image using the dims function from DimensionalData:

dims(slice2)

X Sampled{Int64} 800:1600 ForwardOrdered Regular Points,
Y Sampled{Int64} 1:800 ForwardOrdered Regular Points

Named Dimensions

Each dimension of an AstroImage is named. The automatic dimension names are X, Y, Z, Dim{4}, Dim{5}, and so on; however you can pass in other names or orders to the load function and/or AstroImage contructor:

julia> img = load("eagle-656nmos.fits", 1, (Y,Z))
1600×1600 AstroImage{Float32,2} with dimensions:
  Y Sampled 1:1600 ForwardOrdered Regular Points,
  Z Sampled 1:1600 ForwardOrdered Regular Points

Other useful dimension names are Spec for spectral axes, Pol for polarization data, and Ti for time axes. These are tracked the same way as the automatic dimension names and interact smoothly with any WCS headers. You can give a dimension an arbitrary name using Dim{Symbol}, e.g., Dim{:Velocity}.

You can access AstroImages using dimension names:

eagle[X=100]

1600-element AstroImage{Float32,1} with dimensions: 
  Y Sampled{Int64} 1:1600 ForwardOrdered Regular Points
and reference dimensions: 
  X Sampled{Int64} 100:100 ForwardOrdered Regular Points
    1  0.0
    2  0.0
    3  0.0
    4  0.0
    5  0.0
    ⋮  
 1597  0.0
 1598  0.0
 1599  0.0
 1600  0.0

When indexing into a slice out of a larger parent image or cube, this named access refers to the parent dimensions:

slice1 = eagle[600:800,600:800]
slice1[X=At(700),Y=At(700)] == eagle[X=At(700),Y=At(700)] == eagle[700,700]

true

Cubes

Let's see how this works with a 3D cube.

using AstroImages

HIcube = load(download("http://www.astro.uvic.ca/~wthompson/astroimages/fits/HIdat.fits"))

150×150×450 AstroImage{Float32,3} with dimensions: 
  X Sampled{Int64} Base.OneTo(150) ForwardOrdered Regular Points,
  Y Sampled{Int64} Base.OneTo(150) ForwardOrdered Regular Points,
  Z Sampled{Int64} Base.OneTo(450) ForwardOrdered Regular Points
[:, :, 1]
        1    2    3    4    5    6    7  …  144  145  146  147  148  149  150
   1  NaN  NaN  NaN  NaN  NaN  NaN  NaN     NaN  NaN  NaN  NaN  NaN  NaN  NaN
   2  NaN  NaN  NaN  NaN  NaN  NaN  NaN     NaN  NaN  NaN  NaN  NaN  NaN  NaN
   3  NaN  NaN  NaN  NaN  NaN  NaN  NaN     NaN  NaN  NaN  NaN  NaN  NaN  NaN
   4  NaN  NaN  NaN  NaN  NaN  NaN  NaN     NaN  NaN  NaN  NaN  NaN  NaN  NaN
   ⋮                        ⋮            ⋱    ⋮                        ⋮  
 146  NaN  NaN  NaN  NaN  NaN  NaN  NaN     NaN  NaN  NaN  NaN  NaN  NaN  NaN
 147  NaN  NaN  NaN  NaN  NaN  NaN  NaN     NaN  NaN  NaN  NaN  NaN  NaN  NaN
 148  NaN  NaN  NaN  NaN  NaN  NaN  NaN     NaN  NaN  NaN  NaN  NaN  NaN  NaN
 149  NaN  NaN  NaN  NaN  NaN  NaN  NaN     NaN  NaN  NaN  NaN  NaN  NaN  NaN
 150  NaN  NaN  NaN  NaN  NaN  NaN  NaN  …  NaN  NaN  NaN  NaN  NaN  NaN  NaN
[and 449 more slices...]

Notice how the cube is not displayed automatically. We have to pick a specific slice:

HIcube[Z=228]

Using implot, the world coordinates are displayed automatically:

implot(HIcube[Z=228], cmap=:turbo)

The plot automatically reflects the world coordinates embeded in the file. It displays the x axis in galactic longitude, the y-axis in galactic latitude, and even shows the curved projection from pixel coordinates to galactic coordinates. The title is automatically set to the world coordinate along the Z axis in units of velocity. It also picks up the unit of the data (Kelvins) to display on the colorbar.

If we pick another slice, the title updates accordingly:

implot(HIcube[Z=308], cmap=:turbo)

This works for other slices through the cube as well:

implot(HIcube[Y=45],  cmap=:turbo, aspectratio=0.3)

Custom Dimensions

julia> img = load("img.fits",1,(Y=1:1600,Z=1:1600))
1600×1600 AstroImage{Float32,2} with dimensions:
  Y Sampled 1:1600 ForwardOrdered Regular Points,
  Z Sampled 1:1600 ForwardOrdered Regular Points

Other useful dimension names are Spec for spectral axes, Pol for polarization data, and Ti for time axes. These are tracked the same was as the automatic dimension names and interact smoothly with any WCS headers.

Often times we have images or cubes that we want to index with physical coordinates where setting up a full WCS transform is overkill. In these cases, it's easier to leverage custom dimensions.

For example, one may wish to

julia> img = load("img.fits",1,(X=801:2400,Y=1:2:3200))
1600×1600 AstroImage{Float32,2} with dimensions:
  X Sampled 801:2400 ForwardOrdered Regular Points,
  Y Sampled 1:2:3199 ForwardOrdered Regular Points
...

Unlike OffsetArrays, the usual indexing remains so img[1,1] is still the bottom left of the image; however, data can be looked up according to the offset dimensions using specifiers:

julia> img[X=Near(2000),Y=1..100]
50-element AstroImage{Float32,1} with dimensions:
  Y Sampled 1:2:99 ForwardOrdered Regular Points
and reference dimensions:
  X Sampled 2000:2000 ForwardOrdered Regular Points
  0.0

You can adjust the center of an image's dimensions using recenter:

eagle_cen = recenter(eagle, 801, 801);

Unlike an OffsetArray, eagle_cen[1,1] still refers to the bottom left of the image. This also has no effect on broadcasting; eagle_cen .+ ones(1600,1600) is perfectly valid. However, we see the new centered dimensions when we go to plot the image:

implot(eagle_cen, wcsticks=false)

And we can query positions using the offset dimensions:

implot(eagle_cen[X=-300..300, Y=-300..300], wcsticks=false)