magcutout: Image Cutout Utilities

Description

Functions to subset both raw images and images with associated WCS systems.

Usage

magcutout(image, loc = dim(image)/2, box = c(100, 100), shiftloc = FALSE, paddim = TRUE,
plot = FALSE, ...)
magcutoutWCS(image, header, loc, box = c(100, 100), shiftloc = FALSE, paddim = TRUE,
plot = FALSE, CRVAL1 = 0, CRVAL2 = 0, CRPIX1 = 0, CRPIX2 = 0, CD1_1 = 1, CD1_2 = 0,
CD2_1 = 0, CD2_2 = 1, coord.type = "deg", sep = ":", loc.type=c('coord','coord'), ...)
magWCSradec2xy(RA, Dec, header, CRVAL1 = 0, CRVAL2 = 0, CRPIX1 = 0, CRPIX2 = 0, CD1_1 = 1,
CD1_2 = 0, CD2_1 = 0, CD2_2 = 1, loc.diff = c(0, 0), coord.type = "deg", sep = ":")
magWCSxy2radec(x, y, header, CRVAL1 = 0, CRVAL2 = 0, CRPIX1 = 0, CRPIX2 = 0, CD1_1 = 1,
CD1_2 = 0, CD2_1 = 0, CD2_2 = 1, loc.diff = c(0, 0))

Arguments

image

Numeric matrix; required, the image we want to decorate. If image is a list as created by readFITS, read.fits of magcutoutWCS then the image part of these lists is parsed to image and the correct header part is parsed to header.

header

Full FITS header in table or vector format. Legal table format headers are provided by the read.fitshdr function or the hdr list output of read.fits in the astro package; the hdr output of readFITS in the FITSio package or the header output of magcutoutWCS. If a header is provided then key words will be taken from here as a priority. Missing header keywords are printed out and other header option arguments are used in these cases.

loc

Numeric vector; the [x,y] (magcutout) or [x,y] / [RA,Dec] (magcutoutWCS) location where we want to cutout the image. For magcutoutWCS the unit type can be specified with the loc.type option. Either it is WCS in degrees [RA,Dec] (coord) or pixel [x,y] of the image (image).

box

Numeric vector; the dimensions of the box in pixels (magcutout) or pixels/asec (magcutoutWCS) to cut out from image centred on loc.

Vector or matrix; target right ascension in degrees. If matrix then the first column will be used as RA and the second column as Dec.

Dec

Vector; target declination in degrees. Ignored if RA is a matrix.

Vector or matrix; target x-pixel. If Matrix then the first column will be used as the x-axis and the second column as y-axis. Note this is the R convention of [x,y] (see Notes) not FITS.

Vector; target y-pixel. Ignored if x is a matrix. Note this is the R convention of [x,y] not FITS (see Notes).

loc.diff

The pixel offset to apply. Only relevant if the image being plotted is a cutout from within a FITS legal image.

shiftloc

Logical; should the cutout centre shift from loc away from the image edge if the desired box size extends beyond the edge of the image?

paddim

Logical; should the cutout be padded with image data until it meets the desired box size (if shiftloc is TRUE) or padded with NAs for data outside the image boundary otherwise?

plot

Logical; should a magimage (magcutout) or magimageWCS (magcutoutWCS) plot of the output be generated?

CRVAL1

FITS header CRVAL1 for the Tan Gnomonic projection system. This is the RA in degrees at the location of CRPIX1.

CRVAL2

FITS header CRVAL2 for the Tan Gnomonic projection system. This is the Dec in degrees at the location of CRPIX2.

CRPIX1

FITS header CRPIX1 for the Tan Gnomonic projection system. This is the x pixel value at the location of CRVAL1.

CRPIX2

FITS header CRPIX2 for the Tan Gnomonic projection system. This is the y pixel value at the location of CRVAL2.

CD1_1

FITS header CD1_1 for the Tan Gnomonic projection system. Change in RA-Tan in degrees along x-Axis.

CD1_2

FITS header CD1_2 for the Tan Gnomonic projection system. Change in RA-Tan in degrees along y-Axis.

CD2_1

FITS header CD2_1 for the Tan Gnomonic projection system. Change in Dec-Tan in degrees along x-Axis.

CD2_2

FITS header CD2_2 for the Tan Gnomonic projection system. Change in Dec-Tan in degrees along y-Axis.

coord.type

The units of loc for magcutoutWCS. Allowed options are 'deg' for degress and 'sex' for sexigesimal (i.e. HMS for RA and DMS for Deg).

sep

When coord.type='sex', sep defines the type of separator used for the HMS and DMS strings (i.e. H:M:S and D:M:S would be sep=':', which is the default). See hms2deg and dms2deg for more details.

loc.type

Character vector; specifies what type of location is being provided. The first element specifies the coordinate type for loc and the second element is the coordinate type for box. Either it is WCS in degrees [RA,Dec] / asec ('coord') or pixel [x,y] of the image ('image').

…

Dots are parsed to either magimage (magcutout) or magimageWCS (magcutoutWCS).

Value

A list containing:

image

Numeric matrix; the cutout region of interst centred around loc.

loc

The new loc vector that tranforms the input loc x and y location to the new cutim coordinates. This is in ProFit coordinates, so these values can be used when, e.g., constructing a ProFit modellist structure.

loc.orig

The original location is provided by the input loc.

loc.diff

The x and y offsets of the cutout compared to the original image, where loc + loc.diff = loc.orig exactly.

xsel

Integer vector; the extracted x pixels from the original image that form cutim.

ysel

Integer vector; the extracted y pixels from the original image that form cutim.

loc.WCS

Only output from magcutoutWCS. Numeric vector of the central RA and Dec of the cutout region in degrees.

scale.WCS

Only output from magcutoutWCS. Numeric scalar value of the pixel scale in asec/pixel.

usr.WCS

Only output from magcutoutWCS. Numeric matrix of the RA/Dec extremes of the image in order BL/TL/BR/TR.

approx.map

Only output from magcutoutWCS. A helper function that will usually approximately map RA and Dec numeric vector inputs to a matrix with columns x and y. This is less accurate than magWCSradec2xy, so use that function is the projection is too extreme. Should work for most N-S aligned data pretty well though, and it saves having to correctly track the appropriate header.

header

Only output from magcutoutWCS. The updated header with the correct WCS for the cutout region. Basically this means CRPIX1.new=CRPIX1.org-loc.diff[1] and CRPIX2.new=CRPIX2.org-loc.diff[2].

Details

These functions are on a level trivial, since it is easy to subset matrices and therefore images within R. However these functions track important properties of the subset region that makes it easy to track its location with respect to the original image. Also, they allow direct plotting of the resultant cutout with the most appropriate image functions. In many cases these functions will be used purely for their plotting side effects.

The shiftloc and paddim control the behaviour of the function in the non-trivial case when the desired box size extendeds beyond the edge of the image. If shiftloc is FALSE (the default behaviour), the cutout is guaranteed to be centred on the pixel specified by loc. Then, if paddim is FALSE, the cutout extends only as far as possible until it reaches the edge of the image; otherwise if paddim is TRUE the cutout image is padded with NAs in regions outside the supplied image (the default behaviour). If shiftloc is TRUE, the centre of the cutout will be shifted. In this case, if paddim is FALSE, the cutout will extend at most half of the supplied box size from the given loc; otherwise if paddim is TRUE the cutout will be expanded until it reaches the desired box size or spans the entire image.

Note that if shiftloc is TRUE and paddim is FALSE, the cutout can be larger than box; otherwise, the cutout is guaranteed to be no larger than the specified box size.

Examples

Run this code

# NOT RUN {
temp=matrix(1:121,11)

#The central value is at:

temp[6,6]

print(magcutout(temp, dim(temp)/2, box=c(3,3))$image)

#Given we cutout around the centre of the central pixel [5.5,5.5], the new centre
#relative to the cutout image output should be at [1.5,1.5]:

print(magcutout(temp, dim(temp)/2, box=c(3,3))$loc.orig)
print(magcutout(temp, dim(temp)/2, box=c(3,3))$loc)

# A simple WCS cutout example:

image=readFITS(system.file("extdata", 'VIKING/mystery_VIKING_Z.fits', package="ProFit"))

par(mar=c(3.1,3.1,1.1,1.1))
magcutout(image$imDat, loc=c(50.5,50.5), plot=TRUE)
magcutoutWCS(image, loc=c(50.5,50.5), loc.type = 'image', plot=TRUE)

#We can cut out using the coordinates instead:
print({tempcoord=magWCSxy2radec(50.5,50.5,header=image$hdr)})
magcutoutWCS(image, loc=tempcoord, plot=TRUE)

#You can select coordinates too:

magcutoutWCS(image, loc=c(352.29167, -31.827777), box=c(30,30), plot=TRUE)$loc.WCS
magcutoutWCS(image, loc=c("23:29:10", "-31:49:40"), box=c(30,30) , coord.type = 'sex',
plot=TRUE)$loc.WCS

#We can add RA Dec specific decorations easily:

cutim=magcutoutWCS(image, loc=c(352.2918, -31.82652), box=c(30,30), plot=TRUE)

#Approx overlays:

points(cutim$approx.map(c(352.2918, 352.2897), c(-31.82652, -31.8252)), pch=3, col='blue')

#Exact overlays:

points(magWCSradec2xy(c(352.2918, 352.2897), c(-31.82652, -31.8252), header=cutim$header),
col='red')

#Given we correctly modify the header, we can actually use the cut down image directly:

magimageWCS(cutim)

# Now test the various cutout size options by asking for a large cutout near the image boundary

loc = c(300,340)
box = c(200,200)
loc.type = c("image","image")

# By default, the cutout is exactly the request size, but the centre is shifted:

cutim=magcutoutWCS(image, loc=loc, box=box, coord.type="image", loc.type=loc.type,
	plot=TRUE, shiftloc=TRUE, paddim=TRUE)

# Setting shiftloc=FALSE pads the image with NAs instead, preserving the size:

cutim=magcutoutWCS(image, loc=loc, box=box, coord.type="image", loc.type=loc.type,
	plot=TRUE, shiftloc=FALSE, paddim=TRUE)

# Setting paddim=FALSE returns the largest possible cutout within the image bounds,
# without shifting the centre:

cutim=magcutoutWCS(image, loc=loc, box=box, coord.type="image", loc.type=loc.type, plot=TRUE,
	shiftloc=FALSE, paddim=FALSE)

# Setting paddim=FALSE and shiftloc=TRUE returns a larger cutout, but with at most
# box/2 padding on either side:

cutim=magcutoutWCS(image, loc=loc, box=box, coord.type="image", loc.type=loc.type, plot=TRUE,
	shiftloc=TRUE, paddim=FALSE)

# Setting shiftloc=FALSE and requesting a box size larger than the image returns a cutout with the
# requested box size:

box = c(400,400)
cutim=magcutoutWCS(image, loc=loc, box=box, coord.type="image", loc.type=loc.type, plot=TRUE,
	shiftloc=FALSE, paddim=TRUE)

# }