""" WARNING: If you have ever reason to run this again, first update it to the current keyword set. See the sibling file updateHeaders for what has changed. A preprocessor for the Milcho Tsvetkov Kapteyn field scans. These are bad in a large number of ways; this processor does a first cleanup, resulting in cleaned-up files in data/fits. This file will also insert all non-WCS related header fields. Since the original files are broken as far as pyfits is concerned, I override the process method of the HeaderProcessor. Also, we actually copy the files. This processor uses information prepared by Milcho in res/; this concerns the specification of cutouts, remarks, and lists of observed fields. As we currently ignore rotated images, it's normal that quite a few files remain unprocessed as far as this processor is concerned. Here's positions from the original observation plans to put into RA-ORIG and DEC-ORIG if we feel bored: Kapteyn Special Areas, positions B1900.0 RA.1900.0 D.1900.0 2 19 08 51.3 -00 02.7 3 18 32 24.9 +07 26.8 4 18 56 31.4 +04 07.0 5 19 36 22.4 +10 27.0 6 19 19 37.5 +17 49.4 7 19 33 31.4 +29 23.4 24 06 29 54.9 +10 55.3 25 06 35 50.1 +09 18.4 Kapteyn-Pritchard Area, positions B1900.0 RA.1900.0 D.1900.0 1 00 36 01.9 08 48.6 2 02 50 52.8 07 58.8 3 04 33 41.1 07 40.3 4 08 20 33.0 07 53.4 5 10 55 33.9 06 38.3 6 13 30 53.0 08 48.2 7 15 41 35.6 07 40.0 8 17 58 33.5 08 25.0 9 20 23 15.5 08 06.4 10 22 24 08.5 08 37.1 Kapteyn Selected Area guide stars RA.1900.0 D.1900.0 1900.00 2000.00 68 00 11 00 15 20 00 GSC1179-298 - HD1213 BD+15 30 00 11 25.16 + 15 16 52.9 00 16 34.748 +15 50 13.92 69 01 15 00 15 10 00 GSC1197-472 - HD8110 BD+14 204 01 15 21.45 + 15 10 16.1 01 20 41.101 +15 41 45.74 RA.1900.0 D.1900.0 Epoch 1900.0 Epoch 2000.00 70 02 16 00 15 00 00 GSC1215-1148 02 21 47.135 +15 31 09.09 022147+153109 71 03 11 00 15 00 00 GSC1232-725 03 17 11.624 +15 24 23.86 031712+152424 72 04 10 00 15 10 00 GSC1251-128 04 15 46.284 +15 24 02.50 041546+152402 73 05 14 00 15 00 00 GSC1296-591 05 20 35.750 +15 03 33.53 052036+150334 74 06 15 00 15 10 00 GSC1315-1560 06 20 57.801 +15 03 10.52 062058+150311 75 07 14 00 15 10 00 GSC775-282 07 20 37.877 +14 53 57.95 072038+145358 76 08 15 00 15 00 00 GSC807-359 08 21 09.924 +14 42 27.22 082110+144227 77 09 10 00 14 30 00 BD+14 2057,HD79726 09 16 07.683 +14 07 33.16 091608+140733 78 10 13 00 15 10 00 GSC843-152 10 19 16.756 +14 40 04.91 101917+144005 79 11 17 00 14 50 00 GSC861-423 11 22 54.741 +14 19 22.78 112255+141923 80 12 12 00 15 00 00 GSC879-828 12 17 36.185 +14 26 34.20 121736+142634 81 13 12 00 14 40 00 GSC897-418 13 16 58.485 +14 10 18.31 131658+141018 82 14 14 00 15 20 00 GSC915-14 14 19 26.342 +14 56 08.62 141926+145609 83 15 09 00 14 50 00 BD+1428050,HD135264 15 13 47.155 +14 26 54.32 151347+142654 84 16 12 00 15 00 00 GSC958-289 16 16 46.041 +14 47 25.48 161646+144725 85 17 10 00 15 00 00 GSC990-508 17 15 25.698 +14 54 40.90 171526+145441 86 18 11 00 15 00 00 GSC1568-445 18 16 16.221 +15 04 17.05 181616+150417 87 19 11 00 15 00 00 GSC1599-3918 19 15 48.282 +15 13 20.50 191548+151320 88 20 10 00 15 10 00 BD+14 4240,HD192686 20 15 33.681 +15 26 12.80 201534+152613 89 21 08 00 15 10 00 BD+14 4556,HD202017 21 12 53.927 +15 35 22.69 211254+153523 90 22 12 00 15 10 00 GSC1682-186 22 17 37.348 +15 43 16.65 221737+154317 91 23 13 00 15 00 00 GSC1713-1225 23 18 08.892 +15 36 12.80 231809+153613 """ import array import datetime import math import os import re import sys import ephem from gavo import api from gavo import utils from gavo.grammars import columngrammar from gavo.helpers import processing from gavo.helpers import fitstricks from gavo.utils import fitstools from gavo.utils import pyfits MS = api.makeStruct OBSLOG_GRAMMAR = api.parseFromString(columngrammar.ColumnGrammar, """ coll: 1-6 plate: 8-13 ra: 15-20 dec: 21-27 datetime: 29-42 object: 44-63 exposure: 72-75 emulsion: 76-95 """) def _insertCommentKWs(headerBytes): # insert missing comment keywords where they're necessary result = [] for offset in range(0, len(headerBytes), 80): line = headerBytes[offset:offset+80] if line.startswith("END "): break if line.startswith(" "): line = "COMMENT "+line[8:] result.append(line) return "".join(result)+headerBytes[offset:] def getMinMaxSignedShortSwap(f): """returns the min and max for foreign-endian shorts read from f. """ minVal, maxVal = 1000000, -1000000 while True: buf = f.read(500000) if not buf: break dat = array.array('h') dat.fromstring(buf) dat.byteswap() minVal = min(minVal, min(dat)) maxVal = max(maxVal, max(dat)) return minVal, maxVal def normalizeWhitespace(s): return " ".join(s.split()).strip() # this is what cutout specs look like in our .cutout.txts: CUTOUT_RE = re.compile(r"@\((\d+),(\d+);(\d+),(\d+)\)") # (there can be multiple of these, in which case I'm only using the first one) def loadCutoutMeta(inF, cutouts, comments): """reads cutout meta from inF, dumping the result in cutouts and comments. This comes in a custom text format not worth commenting. """ for ln in inF: if ln.startswith("#") or not ln.strip(): continue series, index, payload = ln.split(" ", 2) mat = CUTOUT_RE.search(payload) if not mat: # if no cutout is given, it's a "plate not found" or similar continue fName = "%s_%s.fits"%(series, index) cutouts[fName] = tuple(int(s) for s in mat.groups()) comments[fName] = normalizeWhitespace( CUTOUT_RE.sub(" ", payload)) def loadObservedMeta(inF, observed): """parses the list of objects observed from inF, dumping the result in observed. """ for ln in inF: if ln.startswith("#") or not ln.strip(): continue series, index, payload = ln.split(" ", 2) fName = "%s_%s.fits"%(series, index) observed[fName] = normalizeWhitespace(payload) def _computeHorizonpars(observer, ra, dec): """returns airmass, zenith distance and hour angle of ra, dec for observer. observer is a pyephem object. """ def sec(x): return 1/math.cos(x) object = ephem.readdb("Star,f,%s,%s,,2000"%( utils.degToHms(ra), utils.degToDms(dec))) object.compute(observer) st = observer.sidereal_time() z = math.pi/2-object.alt airmass = sec(z)*(1-0.0012*(sec(z)**2-1)) return airmass, z/math.pi*180, st*180/math.pi-ra class FixerCopier(processing.HeaderProcessor): def _createAuxiliaries(self, dd): self.destDir = os.path.join(dd.rd.resdir, "data", "fits") if not os.path.isdir(self.destDir): os.mkdir(self.destDir) self.cutouts, self.comments, self.observed = {}, {}, {} with open(os.path.join(dd.rd.resdir, "res", "POT015.cutout.txt")) as f: loadCutoutMeta(f, self.cutouts, self.comments) with open(os.path.join(dd.rd.resdir, "res", "POT080.cutout.txt")) as f: loadCutoutMeta(f, self.cutouts, self.comments) with open(os.path.join(dd.rd.resdir, "res", "POT015.observed.txt")) as f: loadObservedMeta(f, self.observed) with open(os.path.join(dd.rd.resdir, "res", "POT080.observed.txt")) as f: loadObservedMeta(f, self.observed) tmp = list(OBSLOG_GRAMMAR.parse( os.path.join(dd.rd.resdir, "res", "POT015.obslog.txt"))) self.logs = dict( ("%s_%s.fits"%(row["coll"], row["plate"]), row) for row in tmp) tmp = list(OBSLOG_GRAMMAR.parse( os.path.join(dd.rd.resdir, "res", "POT080.obslog.txt"))) self.logs.update(dict( ("%s_%s.fits"%(row["coll"], row["plate"]), row) for row in tmp)) self.obs = ephem.Observer() self.obs.lat, self.obs.long = 52.1934*utils.DEG, 13.0643*utils.DEG self.obs.elevation = 90 def getDestName(self, srcName): return os.path.join(self.destDir, os.path.basename(srcName)) def _isProcessed(self, srcName): return os.path.exists(self.getDestName(srcName)) def _makeNewHeader(self, oldHeader, key): """ """ row = self.logs[key] centerRA = utils.hmsToDeg(row["ra"], sepChar="") centerDec = utils.dmsToDeg(row["dec"], sepChar="") try: startDT = datetime.datetime.strptime(row["datetime"], "%Y%m%d%H%M%S") timeAvailable = True except ValueError: # no time available startDT = datetime.datetime.strptime(row["datetime"], "%Y%m%d") timeAvailable = False if row["exposure"]: exposure = float(row["exposure"])*60 else: exposure = None if timeAvailable: endDT = startDT+datetime.timedelta(seconds=exposure) midDT = startDT+datetime.timedelta(seconds=exposure/2) else: endDT = midDT = None self.obs.date = midDT or startDT airmass, ha, zd = _computeHorizonpars(self.obs, centerRA, centerDec) if "prism" in self.comments[key]: method = "objective prism" prism = "P2" prismang = "7:28" dispers = 350 else: method = "direct photograph" prism = None prismang = None dispers = None if key.startswith("POT015"): telaper = 0.15 scanner = "EPSON Expression 10000XL" pixsize = 10.5833 scanres = 2400 scansoft= 'VueScan, http://www.hamrick.com' wedge = 'TG21S Danes Picta' platesz = 20 fov = 7.65 elif key.startswith("POT080"): telaper = 0.80 scanner = "EPSON Expression V700" pixsize = 15.875 scanres = 1600 scansoft= 'VueScan, http://www.hamrick.com' wedge = 'T13 Danes Picta' platesz = 16 fov = 0.76 else: raise api.ReportableError("Unknown WFPDB id '%s'"%key) res = fitstricks.makeHeaderFromTemplate(fitstricks.WFPDB_TEMPLATE, BITPIX = 16, NAXIS1 = 1, NAXIS2 = 1, OBJECT = row["object"], OBJTYPE = "Kapteyn series", EXPTIME = exposure, NUMEXP = 1, OBSERVAT = "AO Potsdam (42)", SITENAME = "Potsdam Telegrafenberg", SITELONG = self.obs.long/utils.DEG, SITELAT = self.obs.lat/utils.DEG, SITEELEV = self.obs.elevation, TELESCOP = oldHeader["TELESCOP"], DETNAM = "Photographic Plate", TELAPER = oldHeader["TELAPER"], TELFOC = oldHeader["TELFOC"], TELSCALE = oldHeader["TELSCALE"], METHOD = method, FILTER = "none", PRISM = prism, PRISMANG = prismang, DISPERS = dispers, OBSERVER = oldHeader["OBSERVER"], PLATENUM = oldHeader["PLATENUM"], SERIES = key[:6], WFPDB_ID = oldHeader["PLATE-ID"], PLATESZ1 = platesz, PLATESZ2 = platesz, FOV1 = fov, FOV2 = fov, PLATNOTE = self.comments[key], DATE_OBS = utils.formatISODT(startDT), DATE_AVG = utils.formatISODT(midDT) if midDT else None, DATE_END = utils.formatISODT(endDT) if endDT else None, YEAR = api.dateTimeToJYear(startDT), YEAR_AVG = api.dateTimeToJYear(midDT) if midDT else None, JD_AVG = api.dateTimeToJdn(midDT) if midDT else None, RA = utils.degToHms(centerRA), DEC = utils.degToDms(centerDec), RA_DEG = centerRA, DEC_DEG = centerDec, AIRMASS = airmass, ZD = zd, HA = ha, SCANNER = scanner, SCANRES1 = scanres, SCANRES2 = scanres, PIXSIZE1 = pixsize, PIXSIZE2 = pixsize, SCANSOFT= 'VueScan, http://www.hamrick.com', WEDGE = wedge, DATAMIN = 0, DATAMAX = 0, SCANFOC = 'glass', SCANGAM = 1.0, DATESCAN = oldHeader["DATE-SCN"], SCANAUTH = oldHeader["AUTHOR"], ORIGIN = "GAVO Data Center", FILENAME = key, FN_WEDGE = key.replace(".fits", "w.fits"), DATE = utils.formatISODT(datetime.datetime.utcnow())) res.add_comment("This scan originated at Leibniz-Institut" " fuer Astrophysik Potsdam and the Bulgarian Academy of Sciences" " under DFG grant STE 710-6-1 2009-2012") return res def process(self, srcName): if not self.opts.reProcess and self._isProcessed(srcName): return key = os.path.basename(srcName) cuts = self.cutouts.get(key) if cuts is None: return hdu = pyfits.open(srcName)[0] hdu.header = self._makeNewHeader(hdu.header, key) hdu.data = hdu.data[cuts[1]:cuts[3], cuts[0]:cuts[2]] hdu.header.set("DATAMIN", min(hdu.data.ravel())) hdu.header.set("DATAMAX", max(hdu.data.ravel())) # NAXISn updated by pyfits # hdu.scale('int16', '', bzero=32768, bscale=1) hdu.writeto(self.getDestName(srcName), clobber=True) del hdu # broken after scale # The following stuff was orginially done in-place. If we get new issues # of the images, let's hope we won't need it any more # # First fix the broken comment; these are missing COMMENT key words # # after the first one. # with open(srcName, "r+") as f: # f.seek(0) # headerBytes = fitstools.readHeaderBytes(f) # fixedHeaderBytes = _insertCommentKWs(headerBytes) # if len(fixedHeaderBytes)!=len(headerBytes): # raise Exception("Wow, bad bug.") # f.seek(0) # f.write(fixedHeaderBytes) if __name__=="__main__": processing.procmain(FixerCopier, "kapteyn/q", "process_incoming")