This procedure creates a smoothed map of the suface density of
objects on the sky or in a plane.
image processing, catalog processing
DENSITYMAP, x, y, num, map, errmap, head, [cdelt = cdelt, out =
out, /GALACTIC, VERBOSE = verbose, /ROBUST, /DEBUG, /CARTESIAN]
x: Vector of x-coordinate positions for a group of objects, in degrees
y: Corresponding y-coordinate positions for these
objects, in degrees.
num: The number of objects to use when smoothing. See PROCEDURE
section for details.
map: A named variable to hold the density map
errmap: A named variable to hold the error map. This is just map /
sqrt(num - 1)
head: A named variable to hold the fits header describing the map
map, errmap, and head are updated with the correct information.
cdelt: The pixel size of the smoothed map, in degrees /
pixel. Either a scalar or two element array. The default
pixel size is roughly 1/2 the mean input point spacings.
crval: The output map center. The default is the middle of the
input data. A two element array.
naxis: The output map size, in pixels. The default is the smallest
necessary to overlap all input data, given the values for
cdelt and crval. A scalar or two element array.
out: A name used when writing out files. The map and error map
will be written to the file out_map.fits and out_err.fits
galactic: If set, then the input is assumed to be in galactic
verbose: If set, this procedure will print information as it
works. Set to a value 1-4 to control the amount of output
robust: If set, the distances between objects will be computed
using rigorous spherical trigonometry. This will be MUCH slower, as
this appoach is not vectorized well. The default behavior is to
treat the coordinates as cartesian, which is adequate for small
debug: Keyword which runs both robust and non-robust methods, and
reports for how many pixels the non-robust method gets the
cartesian: By default, the procedure assumes that the input
coordinates are on the sky. Set the CARTESIAN keyword to treat the
coordinates as on a flat plane.
The method is taken from section 3.5 of Gutermuth et al
2005ApJ...632..397G. For each pixel, the procedure finds
the distance r to the (num)'th closest object. The surface
density of that pixel is num / (pi r^2).
x = randomn(seed, 500)
y = randomn(seed, 500)
num = 10
DENSITYMAP, x, y, num, map, errmap, head, /verbose, /cartesian
tvscl, congrid(map, 500, 500)
Written by: Chris Beaumont, December 2008
May 11 2009: Added CARTESIAN keyword. cnb.
June 15 2009: Updated things so that VERBIAGE handles text
June 15 2009: Fixed some bugs related to the CARTESIAN
June 15 2009: Enhancements to cdelt. Added crval and naxis
Oct 8 2009: Fixed a bug which led to integer overflow when
manipulating naxis1, naxis2 keywords. cnb