Skip to content

Commit

Permalink
implemented Pluto plates
Browse files Browse the repository at this point in the history
  • Loading branch information
@mommermi
mommermi committed Apr 6, 2018
1 parent 776c9d9 commit cd6c364
Show file tree
Hide file tree
Showing 4 changed files with 1,532 additions and 1,203 deletions.
135 changes: 64 additions & 71 deletions pp_register.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,27 +29,27 @@
import shutil
import logging
import subprocess
import argparse, shlex
import argparse
import shlex
import time
from astropy.io import fits

### pipeline-specific modules
# pipeline-specific modules
import _pp_conf
from catalog import *
import pp_extract
import toolbox
import diagnostics as diag

# only import if Python3 is used
if sys.version_info > (3,0):
if sys.version_info > (3, 0):
from builtins import str

# setup logging
logging.basicConfig(filename = _pp_conf.log_filename,
level = _pp_conf.log_level,
format = _pp_conf.log_formatline,
datefmt = _pp_conf.log_datefmt)

logging.basicConfig(filename=_pp_conf.log_filename,
level=_pp_conf.log_level,
format=_pp_conf.log_formatline,
datefmt=_pp_conf.log_datefmt)


def register(filenames, telescope, sex_snr, source_minarea, aprad,
Expand All @@ -59,9 +59,9 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
registration wrapper
output: diagnostic properties
"""

# start logging
logging.info('starting registration with parameters: %s' % \
logging.info('starting registration with parameters: %s' %
(', '.join([('%s: %s' % (var, str(val))) for
var, val in list(locals().items())])))

Expand All @@ -74,7 +74,7 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
'through pp_prepare?') % filenames[0])
return None

##### run scamp on all image catalogs using different catalogs
# run scamp on all image catalogs using different catalogs
if mancat is not None:
obsparam['astrometry_catalogs'] = [mancat]

Expand All @@ -88,18 +88,18 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,

for cat_idx, refcat in enumerate(obsparam['astrometry_catalogs']):

##### run extract routines
# run extract routines
# ignore saturation: saturated stars are bright and might be necessary
# for SCAMP
if display:
print('* extract sources from %d frames' % len(filenames))

extractparameters = {'sex_snr':sex_snr,
'source_minarea':source_minarea, \
'aprad':aprad, 'telescope':telescope, \
'ignore_saturation':True,
'global_background': True,
'quiet':False}
extractparameters = {'sex_snr': sex_snr,
'source_minarea': source_minarea,
'aprad': aprad, 'telescope': telescope,
'ignore_saturation': True,
'global_background': False,
'quiet': False}

extraction = pp_extract.extract_multiframe(filenames,
extractparameters)
Expand All @@ -110,14 +110,14 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
logging.error('extraction was not successful')
return None

### check if enough sources have been detected in images
# check if enough sources have been detected in images
ldac_files = []
ldac_catalogs = []
for frame in extraction:
if frame['catalog_data'].shape[0] > 10:
ldac_files.append(frame['ldac_filename'])
cat = catalog(frame['ldac_filename'])
cat.read_ldac(frame['ldac_filename'],
cat.read_ldac(frame['ldac_filename'],
frame['fits_filename'],
object_keyword=obsparam['object'],
exptime_keyword=obsparam['exptime'],
Expand All @@ -133,7 +133,7 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
output = {}
fileline = " ".join(ldac_files)

### check if sufficient reference stars are available in refcat
# check if sufficient reference stars are available in refcat
logging.info('check if sufficient reference stars in catalog %s' %
refcat)

Expand All @@ -147,10 +147,11 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,

checkrefcat = catalog(refcat, display=False)
n_sources = checkrefcat.download_catalog(ra, dec,
rad+obsparam['reg_search_radius'],
100, save_catalog=False)
rad +
obsparam['reg_search_radius'],
100, save_catalog=False)
if n_sources < _pp_conf.min_sources_astrometric_catalog:
logging.info(('Only %d sources in astrometric reference catalog; ' \
logging.info(('Only %d sources in astrometric reference catalog; '
+ 'try other catalog') % n_sources)
goodfits, badfits = [], filenames
continue
Expand Down Expand Up @@ -185,8 +186,8 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
'high': '0x00fc'}[source_tolerance]

# assemble arguments for scamp, run it, and wait for it
commandline = 'scamp -c '+obsparam['scamp-config-file']+ \
' -ASTR_FLAGSMASK '+st_code+' -FLAGS_MASK '+st_code+ \
commandline = 'scamp -c '+obsparam['scamp-config-file'] + \
' -ASTR_FLAGSMASK '+st_code+' -FLAGS_MASK '+st_code + \
' -ASTREF_CATALOG FILE' + \
' -ASTREFCAT_NAME ' + refcat + '.cat ' + fileline

Expand All @@ -195,19 +196,19 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
scamp = subprocess.Popen(shlex.split(commandline))
scamp.wait()

##### identify successful and failed WCS registrations based on
##### the contrast values provided by SCAMP
# identify successful and failed WCS registrations based on
# the contrast values provided by SCAMP
scamp = _pp_conf.read_scamp_output()
os.rename('scamp_output.xml', 'astrometry_scamp.xml')
goodfits, badfits = [], []
fitresults = [] # store scamp outputs
fitresults = [] # store scamp outputs
for dat in scamp[1]:
# successful fit
if ((float(dat[scamp[0]['AS_Contrast']]) <
_pp_conf.scamp_as_contrast_limit)
or (float(dat[scamp[0]['XY_Contrast']]) <
_pp_conf.scamp_xy_contrast_limit)
or len(dat) == 0):
or len(dat) == 0):
filename = dat[scamp[0]['Catalog_Name']]
for file in os.listdir('.'):
if file.find(filename[:filename.find('.ldac')]+'.fit') \
Expand All @@ -223,38 +224,37 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
filename = file
goodfits.append(filename)

fitresults.append(\
fitresults.append(
[filename,
float(dat[scamp[0]['AS_Contrast']]),
float(dat[scamp[0]['XY_Contrast']]),
float(dat[scamp[0]['AstromSigma_Reference']].split()[0]),
float(dat[scamp[0]['AstromSigma_Reference']].split()[1]),
float(dat[scamp[0]['Chi2_Reference']]), \
float(dat[scamp[0]['Chi2_Reference']]),
float(dat[scamp[0]['Chi2_Internal']])])

open('registration_succeeded.lst', 'w').writelines("%s\n" %
'\n'.join(goodfits))
if len(goodfits)==0 and len(badfits)==0:
'\n'.join(goodfits))
if len(goodfits) == 0 and len(badfits) == 0:
badfits = filenames
open('registration_failed.lst', 'w').writelines("%s\n" %
'\n'.join(badfits))
'\n'.join(badfits))

output['goodfits'] = goodfits
output['badfits'] = badfits
output['fitresults'] = fitresults
output['catalog'] = refcat
output['goodfits'] = goodfits
output['badfits'] = badfits
output['fitresults'] = fitresults
output['catalog'] = refcat

# check registration outcome

##### check registration outcome

logging.info(' > match succeeded for %d/%d images' % \
logging.info(' > match succeeded for %d/%d images' %
(len(goodfits), len(filenames)))
print('\n################################# ' + \
'REGISTRATION SUMMARY:\n###')
print('### %d/%d images have been registered successfully' % \
(len(goodfits), len(filenames)))
print('###\n###############################' + \
'#######################\n')
print('\n################################# ' +
'REGISTRATION SUMMARY:\n###')
print('### %d/%d images have been registered successfully' %
(len(goodfits), len(filenames)))
print('###\n###############################' +
'#######################\n')

# # registration succeeded for all images
# if len(badfits) == 0:
Expand All @@ -279,22 +279,23 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,

n_success_last_iteration = len(goodfits)

##### update image headers with wcs solutions where registration
##### was successful
# update image headers with wcs solutions where registration
# was successful
logging.info('update image headers with WCS solutions ')

for filename in goodfits:
# remove fake wcs header keys
fake_wcs_keys = ['RADECSYS', 'CTYPE1', 'CTYPE2', 'CRVAL1', 'CRVAL2',
'CRPIX1', 'CRPIX2', 'CD1_1', 'CD1_2', 'CD2_1',
'CD2_2', 'RADESYS']
fake_wcs_keys = ['RADECSYS', 'CTYPE1', 'CTYPE2', 'CRVAL1', 'CRVAL2',
'CRPIX1', 'CRPIX2', 'CD1_1', 'CD1_2', 'CD2_1',
'CD2_2', 'RADESYS']
hdu = fits.open(filename, mode='update', verify='silentfix',
ignore_missing_end=True)
for fake_key in fake_wcs_keys:
hdu[0].header[fake_key] = ''

# read new header files
newhead = open(filename[:filename.find('.fit')]+'.head','r').readlines()
newhead = open(filename[:filename.find(
'.fit')]+'.head', 'r').readlines()

for line in newhead:
key = line[:8].strip()
Expand All @@ -305,7 +306,7 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
comment = line[30:].strip()
if key.find('END') > -1:
break
#print key, '|', value, '|', comment
# print key, '|', value, '|', comment
hdu[0].header[key] = (str(value), comment)

# other header keywords
Expand All @@ -320,31 +321,29 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
if len(goodfits) == len(filenames):
os.remove(filename[:filename.find('.fit')]+'.head')


if len(badfits) == len(filenames):
if display:
print('ERROR: registration failed for all images')
logging.error('ERROR: registration failed for all images')
return output

##### print astrometry output file
# print astrometry output file
outf = open('best_astrometry.dat', 'w')
outf.writelines('# filename AS_contrast XY_contrast ' \
outf.writelines('# filename AS_contrast XY_contrast '
+ 'Chi2_catalog Chi2_int Pos_uncertainty(arcsec)\n')
for idx, data in enumerate(fitresults):
outf.writelines('%25.25s %5.2f %5.2f %10.7f %10.7f %7.4f\n' % \
outf.writelines('%25.25s %5.2f %5.2f %10.7f %10.7f %7.4f\n' %
(data[0], data[1], data[2], data[5], data[6],
numpy.sqrt(data[3]**2+data[4]**2)))
outf.close()


### extraction output
# extraction output
#
# -> see pp_extract.py
#
##

### output content
# output content
#
# { 'good_fits' : list of fits where registration succeeded,
# 'bad fits' : list of fits where registration failed,
Expand All @@ -353,17 +352,15 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
# }
###

##### create diagnostics
# create diagnostics
if diagnostics:
diag.add_registration(output, extraction)
diag.add_registration(output, extraction)

logging.info('Done! -----------------------------------------------------')

return output




if __name__ == '__main__':

# command line arguments
Expand All @@ -386,8 +383,7 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
source_tolerance = args.source_tolerance
filenames = args.images


### read telescope and filter information from fits headers
# read telescope and filter information from fits headers
# check that they are the same for all images
instruments = []
for filename in filenames:
Expand All @@ -402,7 +398,6 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
raise KeyError('cannot identify telescope/instrument; please update'
'_pp_conf.instrument_keys accordingly')


# assign telescope parameters (telescopes.py)
telescope = _pp_conf.instrument_identifiers[instruments[0]]
obsparam = _pp_conf.telescope_parameters[telescope]
Expand All @@ -425,5 +420,3 @@ def register(filenames, telescope, sex_snr, source_minarea, aprad,
source_minarea, aprad, mancat, obsparam,
source_tolerance,
display=True, diagnostics=True)


0 comments on commit cd6c364

Please sign in to comment.