In [1]:
!ls *fits
wdd7.040920_0452.051_6.fits wdd7.080104_0214.1025_6.fits
In [2]:
import astropy.io.fits as afits
from astropy.wcs import WCS
from astropy.visualization import ZScaleInterval
import matplotlib
%matplotlib notebook
%pylab
Using matplotlib backend: nbAgg
Populating the interactive namespace from numpy and matplotlib
In [3]:
f1 = afits.open('wdd7.040920_0452.051_6.fits')
f2 = afits.open('wdd7.080104_0214.1025_6.fits')
f1wcs = WCS(f1[0].header)
f2wcs = WCS(f2[0].header)
zscaler = ZScaleInterval(nsamples=1000, contrast=0.25)
f1d = zscaler(f1[0].data)
f2d = zscaler(f2[0].data)
fig = figure(figsize=(10,4))
ax = fig.add_subplot(1,1,1)
ax.imshow(f1d.T, cmap='Reds')
ax.imshow(f2d.T, cmap='Blues', alpha=0.5)
tight_layout()
xlabel('x')
ylabel('y')
savefig('out/original_misalignment.pdf')
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Default coordinate system
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
In [4]:
!solve-field -h
This program is part of the Astrometry.net suite.
For details, visit http://astrometry.net.
Git URL https://github.com/dstndstn/astrometry.net
Revision 0.73, date Thu_Nov_16_08:30:44_2017_-0500.
Usage: solve-field [options] [<image-file-1> <image-file-2> ...] [<xyls-file-1> <xyls-file-2> ...]
You can specify http:// or ftp:// URLs instead of filenames. The "wget" or "curl" program will be used to retrieve the URL.
Options include:
-h / --help: print this help message
-v / --verbose: be more chatty -- repeat for even more verboseness
-D / --dir <directory>: place all output files in the specified directory
-o / --out <base-filename>: name the output files with this base name
-b / --backend-config <filename>: use this config file for the
"astrometry-engine" program
--config <filename>: use this config file for the "astrometry-engine" program
--batch: run astrometry-engine once, rather than once per input file
-f / --files-on-stdin: read filenames to solve on stdin, one per line
-p / --no-plots: don't create any plots of the results
--plot-scale <scale>: scale the plots by this factor (eg, 0.25)
--plot-bg <filename (JPEG)>: set the background image to use for plots
-G / --use-wget: use wget instead of curl
-O / --overwrite: overwrite output files if they already exist
-K / --continue: don't overwrite output files if they already exist; continue
a previous run
-J / --skip-solved: skip input files for which the 'solved' output file
already exists; NOTE: this assumes single-field input files
--fits-image: assume the input files are FITS images
-N / --new-fits <filename>: output filename of the new FITS file containing
the WCS header; "none" to not create this file
-Z / --kmz <filename>: create KMZ file for Google Sky. (requires wcs2kml)
-i / --scamp <filename>: create image object catalog for SCAMP
-n / --scamp-config <filename>: create SCAMP config file snippet
-U / --index-xyls <filename>: output filename for xylist containing the image
coordinate of stars from the index
--just-augment: just write the augmented xylist files; don't run
astrometry-engine.
--axy <filename>: output filename for augment xy list (axy)
--temp-axy: write 'augmented xy list' (axy) file to a temp file
--timestamp: add timestamps to log messages
-7 / --no-delete-temp: don't delete temp files (for debugging)
-L / --scale-low <scale>: lower bound of image scale estimate
-H / --scale-high <scale>: upper bound of image scale estimate
-u / --scale-units <units>: in what units are the lower and upper bounds?
choices: "degwidth", "degw", "dw" : width of the image, in degrees (default)
"arcminwidth", "amw", "aw" : width of the image, in arcminutes
"arcsecperpix", "app": arcseconds per pixel
"focalmm": 35-mm (width-based) equivalent focal length
-8 / --parity <pos/neg>: only check for matches with positive/negative parity
(default: try both)
-c / --code-tolerance <distance>: matching distance for quads (default: 0.01)
-E / --pixel-error <pixels>: for verification, size of pixel positional error
(default: 1)
-q / --quad-size-min <fraction>: minimum size of quads to try, as a fraction
of the smaller image dimension, default: 0.1
-Q / --quad-size-max <fraction>: maximum size of quads to try, as a fraction
of the image hypotenuse, default 1.0
--odds-to-tune-up <odds>: odds ratio at which to try tuning up a match that
isn't good enough to solve (default: 1e6)
--odds-to-solve <odds>: odds ratio at which to consider a field solved
(default: 1e9)
--odds-to-reject <odds>: odds ratio at which to reject a hypothesis (default:
1e-100)
--odds-to-stop-looking <odds>: odds ratio at which to stop adding stars when
evaluating a hypothesis (default: HUGE_VAL)
--use-sextractor: use SExtractor rather than built-in image2xy to find sources
--sextractor-config <filename>: use the given SExtractor config file. Note
that CATALOG_NAME and CATALOG_TYPE values will be over-ridden by
command-line values. This option implies --use-sextractor.
--sextractor-path <filename>: use the given path to the SExtractor executable.
Default: just 'sex', assumed to be in your PATH. Note that you can
give command-line args here too (but put them in quotes), eg:
--sextractor-path 'sex -DETECT_TYPE CCD'. This option implies
--use-sextractor.
-3 / --ra <degrees or hh:mm:ss>: only search in indexes within 'radius' of the
field center given by 'ra' and 'dec'
-4 / --dec <degrees or [+-]dd:mm:ss>: only search in indexes within 'radius'
of the field center given by 'ra' and 'dec'
-5 / --radius <degrees>: only search in indexes within 'radius' of the field
center given by ('ra', 'dec')
-d / --depth <number or range>: number of field objects to look at, or range
of numbers; 1 is the brightest star, so "-d 10" or "-d 1-10" mean look
at the top ten brightest stars only.
--objs <int>: cut the source list to have this many items (after sorting, if
applicable).
-l / --cpulimit <seconds>: give up solving after the specified number of
seconds of CPU time
-r / --resort: sort the star brightnesses by background-subtracted flux; the
default is to sort using acompromise between background-subtracted and
non-background-subtracted flux
-6 / --extension <int>: FITS extension to read image from.
--invert: invert the image (for black-on-white images)
-z / --downsample <int>: downsample the image by factor <int> before running
source extraction
--no-background-subtraction: don't try to estimate a smoothly-varying sky
background during source extraction.
--sigma <float>: set the noise level in the image
--nsigma <float>: number of sigma for a source detection; default 8
-9 / --no-remove-lines: don't remove horizontal and vertical overdensities of
sources.
--uniformize <int>: select sources uniformly using roughly this many boxes
(0=disable; default 10)
--no-verify-uniformize: don't uniformize the field stars during verification
--no-verify-dedup: don't deduplicate the field stars during verification
-C / --cancel <filename>: filename whose creation signals the process to stop
-S / --solved <filename>: output file to mark that the solver succeeded
-I / --solved-in <filename>: input filename for solved file
-M / --match <filename>: output filename for match file
-R / --rdls <filename>: output filename for RDLS file
--sort-rdls <column>: sort the RDLS file by this column; default is ascending;
use "-column" to sort "column" in descending order instead.
--tag <column>: grab tag-along column from index into RDLS file
--tag-all: grab all tag-along columns from index into RDLS file
-j / --scamp-ref <filename>: output filename for SCAMP reference catalog
-B / --corr <filename>: output filename for correspondences
-W / --wcs <filename>: output filename for WCS file
-P / --pnm <filename>: save the PNM file as <filename>
-k / --keep-xylist <filename>: save the (unaugmented) xylist to <filename>
-A / --dont-augment: quit after writing the unaugmented xylist
-V / --verify <filename>: try to verify an existing WCS file
--verify-ext <extension>: HDU from which to read WCS to verify; set this
BEFORE --verify
-y / --no-verify: ignore existing WCS headers in FITS input images
-g / --guess-scale: try to guess the image scale from the FITS headers
--crpix-center: set the WCS reference point to the image center
--crpix-x <pix>: set the WCS reference point to the given position
--crpix-y <pix>: set the WCS reference point to the given position
-T / --no-tweak: don't fine-tune WCS by computing a SIP polynomial
-m / --temp-dir <dir>: where to put temp files, default /tmp
The following options are valid for xylist inputs only:
-F / --fields <number or range>: the FITS extension(s) to solve, inclusive
-w / --width <pixels>: specify the field width
-e / --height <pixels>: specify the field height
-X / --x-column <column-name>: the FITS column containing the X coordinate of
the sources
-Y / --y-column <column-name>: the FITS column containing the Y coordinate of
the sources
-s / --sort-column <column-name>: the FITS column that should be used to sort
the sources
-a / --sort-ascending: sort in ascending order (smallest first); default is
descending order
Note that most output files can be disabled by setting the filename to "none".
(If you have a sick sense of humour and you really want to name your output
file "none", you can use "./none" instead.)
--ra, --dec and --radius: Restrict the solution to some radius around RA and Dec. The regular telescope WCS should be plenty for an initial guess.--scale-units, --scale-low and --scale-high: You might not know the exact pixel scale (and it's a function of where you are on the detector in any case, but you also set limits from this based on the existing CD1_2, CD2_1-D, -N: Write to an output directory and write out a new FITS file with the solved WCS.--parity: You can usually set this and get a speedup of 2xastropy, or just use imhead from the WCSTools package at the command line
In [5]:
print(WCS(f1[0].header))
WCS Keywords
Number of WCS axes: 2
CTYPE : 'RA---TAN' 'DEC--TAN'
CRVAL : 38.355765019813 -8.9280976775107
CRPIX : -1029.2364 4097.4406
CD1_1 CD1_2 : -1.1704629057241e-07 7.40877225586965e-05
CD2_1 CD2_2 : 7.43153648473077e-05 7.04796871082368e-07
NAXIS : 1024 4096
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Default coordinate system
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
In [6]:
!imhead wdd7.040920_0452.051_6.fits
SIMPLE = T / Fits standard
BITPIX = 16 /
NAXIS = 2 / Number of axes
NAXIS1 = 1024 /
NAXIS2 = 4096 /
EXTEND = F / File may contain extensions
BSCALE = 1.000000 /
BZERO = 32768.000000 /
ORIGIN = 'NOAO-IRAF FITS Image Kernel July 2003' / FITS file originator
DATE = '2005-09-04T19:41:28' / Date FITS file was generated
IRAF-TLM= '15:41:27 (04/09/2005)' / Time of last modification
OBJECT = 'wdd7 ' / Name of the object observed
NEXTEND = 1 /
FILENAME= 'wdd7.040920_0452.051' / Original host filename
OBSTYPE = 'object ' / Observation type
PREFLASH= 0.000000 / Preflash time (sec)
RADECSYS= 'FK5 ' / Default coordinate system
RADECEQ = 2000. / Default equinox
RA = '02:33:25.68' / RA of observation (hr)
DEC = '-08:55:29.70' / DEC of observation (deg)
COMMENT 1 blank line
TIMESYS = 'UTC approximate' / Time system
DATE-OBS= '2004-09-20T04:54:34.0' / Date of observation start (UTC approximate)
TIME-OBS= '04:54:34.003' / Time of observation start
MJD-OBS = 53268.20456022 / MJD of observation start
MJDHDR = 53268.20447685 / MJD of header creation
LSTHDR = '00:08:55.7 ' / LST of header creation
COMMENT 1 blank line
OBSERVAT= 'CTIO ' / Observatory
TELESCOP= 'CTIO 4.0 meter telescope' / Telescope
TELRADEC= 'FK5 ' / Telescope coordinate system
TELEQUIN= 2000.0 / Equinox of tel coords
TELRA = '02:33:25.69 ' / RA of telescope (hr)
TELDEC = '-08:55:29.7 ' / DEC of telescope (deg)
HA = '-02:24:44.3 ' / hour angle (H:M:S)
ZD = '39.9 ' / Zenith distance
AIRMASS = 1.303 / Airmass
TELFOCUS= '16914 ' / Telescope focus
CORRCTOR= 'Blanco Corrector' / Corrector Identification
ADC = 'Blanco ADC' / ADC Identification
COMMENT 1 blank line
DETECTOR= 'Mosaic2 ' / Detector
DETSIZE = '[1:8192,1:8192]' / Detector size for DETSEC
NCCDS = 8 / Number of CCDs
NAMPS = 16 / Number of Amplifiers
PIXSIZE1= 15. / Pixel size for axis 1 (microns)
PIXSIZE2= 15. / Pixel size for axis 2 (microns)
PIXSCAL1= 0.270 / Pixel scale for axis 1 (arcsec/pixel)
PIXSCAL2= 0.270 / Pixel scale for axis 2 (arcsec/pixel)
RAPANGL = 0. / Position angle of RA axis (deg)
DECPANGL= 90. / Position angle of DEC axis (deg)
FILPOS = 5 / Filter position
FILTER = 'R Harris c6004 ' / Filter name(s)
SHUTSTAT= 'guide ' / Shutter status
TV1FOC = -1.8930000000000E+00 / North TV Camera focus position
TV2FOC = -1.5000000000000E+00 / South Camera focus position
ENVTEM = 1.3600000000000E+01 / Ambient temperature (C)
DEWAR = 'Mosaic2 Dewar' / Dewar identification
DEWTEM = 0.0000000000000E+00 / Dewar temperature (C)
DEWTEM2 = -6.1400002000000E+01 / Fill Neck temperature (C)
DEWTEM3 = 'N2 88.7'
CCDTEM = -9.5000000000000E+01 / CCD temperature (C)
CCDTEM2 = 'CCD 169.7'
COMMENT 1 blank line
WEATDATE= 'Sep 20 04:51:07 2004' / Date and time of last update
WINDSPD = '16.0 ' / Wind speed (mph)
WINDDIR = '354 ' / Wind direction (degrees)
AMBTEMP = '12.4 ' / Ambient temperature (degrees C)
HUMIDITY= '27 ' / Ambient relative humidity (percent)
PRESSURE= '781 ' / Barometric pressure (millibars)
DIMMSEE = 'mysql_query():Unknown Error seeing=column' / Tololo DIMM seeing
COMMENT 1 blank line
CONTROLR= 'Mosaic Arcon' / Controller identification
CONSWV = '13July99ver7_30' / Controller software version
AMPINTEG= 3600 / (ns) Double Correlated Sample time
READTIME= 14400 / (ns) unbinned pixel read time
ARCONWD = 'Obs Sun Sep 19 15:55:48 2004' / Date waveforms last compiled
COMMENT 1 blank line
OBSERVER= 'armin rest' / Observer(s)
PROPOSER= '<unknown>' / Proposer(s)
PROPOSAL= '<unknown>' / Proposal title
PROPID = '2002B-0007'
OBSID = 'ct4m.20040920T045434' / Observation ID
COMMENT 1 blank line
IMAGESWV= 'mosdca (Jun99), mosaicsrc.tcl (Nov02)' / Image creation software vers
KWDICT = 'MosaicV1.dic (Sep97)' / Keyword dictionary
COMMENT 1 blank line
OTFDIR = 'mscdb$noao/ctio/4meter/caldir/Mosaic2/' / OTF calibration directory
COMMENT 1 blank line
CHECKSUM= '<unknown>' / Header checksum
DATASUM = '<unknown>' / Data checksum
CHECKVER= '<unknown>' / Checksum version
RECNO = 0 / NOAO archive sequence number
ARCONGI = 2 / Gain selection (index into Gain Table)
COMMENT 1 blank line
FIXIMFIL= '2004-05-11.fixim'
WCSDB = '2004-09-08.sm64.sm040915.IP.050213.1042.refcat.db'
WCSCAT = 'W2 '
EWCSXAS = 4.900000000000E-02
EWCSYAS = 4.400000000000E-02
WCSEXE = 'msccmatch'
WCSFLAG = 1.000000000000E+00
WCSNOBJ = 4.080000000000E+02
XTALKFIL= '/pool/cluster/SMSN/mscv6.0/mscpipe/config/CTIO4m/SMSN/xtalk/2004-05-'
IMAGEID = 6 / Image identification
EXPTIME = 200.000 / Exposure time in secs
DARKTIME= 209.076 / Total elapsed time in secs
COMMENT 1 blank line
CCDNAME = 'SITe #98261FACR06-02 (NOAO 28)' / CCD name
AMPNAME = 'SITe #98261FACR06-02 (NOAO 28), lower right (Amp12)' / Amplifier name
GAIN = 2.218672 /
RDNOISE = 11.098980 /
SATURATE= 36100.000000 /
CONHWV = 'ACEB002_AMP12' / Controller hardware version
ARCONG = 1.9 / gain expected for amp 212 (e-/ADU)
BPM = 'maskdir$wdd7.040920_0452.051_6.pl' / Bad pixel mask
CCDSIZE = '[1:2048,1:4096]' / CCD size
CCDSUM = '1 1 ' / CCD pixel summing
CCDSEC = '[1025:2048,1:4096]' / CCD section
AMPSEC = '[1048:25,1:4096]' / Amplifier section
DETSEC = '[5121:6144,1:4096]' / Detector section
COMMENT 1 blank line
ATM1_1 = -1. / CCD to amplifier transformation
ATM2_2 = 1. / CCD to amplifier transformation
ATV1 = 2073.0 / CCD to amplifier transformation
ATV2 = 0. / CCD to amplifier transformation
LTM1_1 = 1.0 / CCD to image transformation
LTM2_2 = 1.0 / CCD to image transformation
LTV1 = -1024.000000 /
DTM1_1 = 1. / CCD to detector transformation
DTM2_2 = 1. / CCD to detector transformation
DTV1 = 4096. / CCD to detector transformation
DTV2 = 0. / CCD to detector transformation
COMMENT 1 blank line
WCSASTRM= 'ct4m.20020927T000352 (USNO N R Harris c6004) by F. Valdes 2002-10-24'
EQUINOX = 2000. / Equinox of WCS
WCSDIM = 2 / WCS dimensionality
CTYPE1 = 'RA---TNX' / Coordinate type
CTYPE2 = 'DEC--TNX' / Coordinate type
CRVAL1 = 38.355765019813 / Coordinate reference value
CRVAL2 = -8.9280976775107 / Coordinate reference value
CRPIX1 = -1029.236400 /
CRPIX2 = 4097.440600 /
CD1_1 = -1.1704629057241E-7 / Coordinate matrix
CD2_1 = 7.43153648473077E-5 / Coordinate matrix
CD1_2 = 7.40877225586965E-5 / Coordinate matrix
CD2_2 = 7.04796871082368E-7 / Coordinate matrix
WAT0_001= 'system=image' / Coordinate system
WAT1_001= 'wtype=tnx axtype=ra lngcor = "3. 4. 4. 2. -0.3037390893371674 -2.198'
WAT1_002= '270304793818E-4 0.06891904927838492 0.154369562810826 0.001032235381'
WAT1_003= '594434 0.01073692175025741 -0.001523675306727922 -0.115624212670156 '
WAT1_004= '-0.002962541566716242 0.008171056275216731 0.01825099713427637 -0.00'
WAT1_005= '408958964870908 -0.1345308755398008 0. "'
WAT2_001= 'wtype=tnx axtype=dec latcor = "3. 4. 4. 2. -0.3037390893371674 -2.19'
WAT2_002= '8270304793818E-4 0.06891904927838492 0.154369562810826 -9.4230844283'
WAT2_003= '14165E-4 -0.003040895251874454 0.002405747773738591 0.00617280507746'
WAT2_004= '1138 0.01021931138769624 -0.01297144901541361 -0.1173773859380067 -0'
WAT2_005= '.02981745508215422 0.05498004218705955 0. "'
XTALKCOR= 'Sep 4 5:35 Crosstalk is 1.66E-4*im7+0.00173*im8'
MSCCMATC= 'Sep 4 13:03 0.00/0.01 1.000/1.000 -0.000/-0.000'
BOOTSTRP= 1 /
BTSTRSCL= 1.126800 /
BTSTRMET= 'extcat'
BTSTRRMS= 0.001100 /
BTSTRN = 10 /
BTSTRFIL= 'sm040909.0x0104.sm041009.bootstrapflats'
FLATNAME= 'frames/0x060104.flat.il.sm040920_6.dome.fits'
BIASNAME= 'frames/0x060100.bias.sm040920_6.fits'
BADPIXNA= 'frames/0x060100.bpm.sm040511_6.fits'
FRINGEFR= 'NONE'
CALDATE = 'Tue_Nov_24_05:58:41_2009'
FLATTYPE= 20 /
SKYADU = 1261.448081 /
PROJNAME= 'ESSENCE/w project'
PROPTITL= 'The w Project: Measuring the Equation of State of the Universe'
PROJPI = 'Suntzeff, N.'
FILEVER = 1.000000 /
PIPEVER = 10.000000 /
SUBDIR = 'sm040920/6' / pipeline subdir of the flattened file
FITSNAME= 'wdd7.040920_0452.051_6.fits' / fits image name
NOISEIM = 0 / flag: set to 1 if noise image exists
MASKIM = 0 / flag: set to 1 if mask image exists
PHOTCODE= '0x60104 '
END
In [7]:
!wcs-to-tan -h
This program is part of the Astrometry.net suite.
For details, visit http://astrometry.net.
Git URL https://github.com/dstndstn/astrometry.net
Revision 0.73, date Thu_Nov_16_08:30:44_2017_-0500.
Usage: wcs-to-tan
-w <WCS input file>
[-e <extension>] FITS HDU number to read WCS from (default 0 = primary)
[-t]: just use TAN projection, even if SIP extension exists.
[-L]: force WCSlib
[-x x-lo]
[-y y-lo]
[-W x-hi]
[-H y-hi]
[-N grid-n]
-o <WCS output file>
[-v]: verbose
In [8]:
!cat Register_images_solutions_commands.txt
solve-field --ra 38.35576501981 --dec -8.9280976775107 --radius 0.5 --scale-units arcsecperpix --scale-low 0.24 --scale-high 0.28 --fits-image wdd7.040920_0452.051_6.fits -N wdd7.040920_0452.051_6.solved.fits --parity pos
wcs-to-tan -w wdd7.040920_0452.051_6.wcs -t -o wdd7.040920_0452.051_6.wcs.tan
solve-field --ra 38.35576501981 --dec -8.9280976775107 --radius 0.5 --scale-units arcsecperpix --scale-low 0.24 --scale-high 0.28 --fits-image wdd7.080104_0214.1025_6.fits -N wdd7.080104_0214.1025_6.solved.fits --parity pos
wcs-to-tan -w wdd7.080104_0214.1025_6.wcs -t -o wdd7.080104_0214.1025_6.wcs.tan
Content source: LSSTC-DSFP/LSSTC-DSFP-Sessions
Similar notebooks: