IO support for reading/writing QDP files

One thing that has been bugging me for quite a long time is that astropy does not support qdp files which are quite common in the X-ray community. It is not difficult to write a small procedure to read them in, but it is just one of those things that keep coming back and while I can read nearly all files using astropy, I can't read the qdp files. For that I need to import a different routine. And yet, these are very commonly used in the X-ray community. I would like to suggest to add a reader/ perhaps writer for those kind of files as well.

Can you point to a spec of the file format?

Yes, here is a reference
https://heasarc.gsfc.nasa.gov/ftools/others/qdp/node4.html

Do you have any interest in writing the small procedure to contribute to astropy?

no, not really, since I would have to learn about astropy internals. Its also a time issue. It should be rather straightforward though, especially for someone who knows astropy.io.

Here is an example program I use to rebin the data in a qdp file. Its main difference to a normal data file is that it can contain multiple tables separated by a 'no no no no no ' line. So reading in a qdp file would give a collection of tables.

def qdp_bin(qdpfile,weights='yes',rate='yes',maxtimebin=0.3,
   minsigma=2.5,maxrate=50.0,multiply=None,chatter=0):
   ''' read qdp file bin up the values using (weighted RMS values)
   
       Parameters
       ==========
       qdpfile : path
          qdp input file with !FILTER records in each section
          where the format is:
          !FILTER filename photometric_system_name
       weights: ["yes","no"]
          average using errors as weights?
       rate : ["yes","no"]
          is the input rates or magnitudes?
       maxtimebin : float
          the maximum size for the time bin you allow
       minsigma : float
          the minimum sigma that is required for a good
          data point.  
       maxrate : float
          the maximumrate that is allowed
       multiply : float
          multiply the magnitude/rate by this factor
       chatter : int
       
   
       Note
       ====
       coadding is based on least squares fit. High rates for uvot 
       need binning using raw source and background data with correction 
       to get the correct error and flux. High bins are flagged.  
       Two-sided errors are converted to one-sided. 
       multiply can be provided as a dictionary for all filters in the qdp file
   '''
   from numpy import log10, sqrt
   #check number of section
   try:
     f = open(qdpfile)
     recs = f.readlines()
     f.close()
   except:
     print 'qdp_bin: input filename error'
     return
   file1 = qdpfile.split('.')[0:-1][0] 
   #    
   # open output file
   #  
   f = open(file1+'_bin.qdp','w')
   #  
   #  set values for bin 
   #
   tmin = None
   tmax = 0
   dt = 0
   x = 0
   dx = 1
   var = 0
   band = ""
   mult = 1.0
   for rec in recs:
      s = rec.split()
      if len(rec) == 0:                    # blank line
         f.write(rec)
      elif rec.split()[0].upper() == '!FILTER':    # define new filter band
         band = rec.split()[1]
         psys = rec.split()[2]
         f.write(rec)
         if multiply != None:
            if multiply[band]:
               mult = multiply[band]
               if chatter > 0: print 'multiplying '+band+' with ',mult
            else:
               mult = 1.0
      elif rec.split()[0][0] == '!':
         f.write(rec)
      elif len(s) > 3:
         #if ((s[0] == 'NO') & (s[1] == 'NO') & (s[2] == 'NO') & (s[3] == 'NO') ):
         if (s[0] == 'NO'):
            # if there is a stored record, output and clear buffer 
            #f.write(str(t)+'  '+str(dt)+' '+str(rx)+' '+str(drx)+' ! '+band+' last point this filter \n' )
            f.write( ("%12.2f  %13.2f  %11.4e %11.4e  ! %-15s  -last point this filter\n" % (t,dt,rx,drx,band)))
            f.write(rec)
            tmin = None ; tmax = 0 ; x = 0 ; dx = 1 ; var = 0    
         else: 
            try: 
               time = float(rec.split()[0])
               halfdt = float(rec.split()[1])   
               rate = float(rec.split()[2])
               rate *= mult
               rerr1 = float(rec.split()[3])
               rerr1 *= mult
               if rate <= 0.:
                  rate = 0.
                  rerr1 = 0. 
                  print 'ignoring '+rec
                  # do not add in
               if tmin == None: tmin = time-halfdt
               tmax = time+halfdt
               t  = 0.5*(tmin+tmax)          # mid time
               dt = 0.5*(tmax-tmin)          # halft the time bin width
               x += rate*2*halfdt            # total #counts in bin   
               
               try:
                  rerr2 = rec.split()[4]
                  if rerr2[0] == '!': 
                     onesided = True
                     var += (2*rerr1*halfdt)**2  # total variance in #counts    ??? wrong I think ???     
                  else:
                     rerr2 = float(rerr2)
                     rerr2 *= mult
                     onesided = False
                     var += (rerr1**2+rerr2**2)*halfdt**2 # total variance ??? wrong I think ??? in #counts
               except RuntimeError:
                  print 'runtime error in qdp_bin'
                  raise
               dx = sqrt(var)                 # error estimate
               sigma = x/dx                   # proxy for standard deviation  
               rx = x/(dt*2.)
               drx = dx/(dt*2.)  
               print "%f, %f, %f - %f,%f - %f, %f, %f, %f " %(tmin,t,tmax, halfdt,dt, rate,rerr1,rx,drx)               
               #
               #  test if the binned values should be written
               #
               if rate > maxrate: 
                  f.write( ("%12.2f  %13.2f  %11.4e %11.4e  ! %-15s WARNING: use correct coi-form when coadding\n" % (t,dt,rx,drx,band)))
                  #f.write(str(t)+'  '+str(dt)+' '+str(rx)+' '+str(drx)+\
                  #                          ' !coadded WARNING: should use correct coi-corrections \n' )
                  tmin = None ; tmax = 0 ; x = 0 ; dx = 1 ; var = 0      
               elif sigma > minsigma: 
                  #f.write(str(t)+'  '+str(dt)+' '+str(rx)+' '+str(drx)+' ! '+band+' sigma > '+str(minsigma)+' \n' )
                  f.write( ("%12.2f  %13.2f  %11.4e %11.4e  ! %-15s  sigma > %-10.3f\n" % (t,dt,rx,drx,band,minsigma)))
                  tmin = None ; tmax = 0 ; x = 0 ; dx = 1 ; var = 0      
               elif dt*2. > maxtimebin*t:       
                  #f.write(str(t)+'  '+str(dt)+' '+str(rx)+' '+str(drx)+' ! '+band+' delta_time > '+str(maxtimebin)+'*(time) \n' )
                  f.write( ("%12.2f  %13.2f  %11.4e %11.4e  ! %-15s  delta time > %-10.3f\n" % (t,dt,rx,drx,band,maxtimebin)))
                  tmin = None ; tmax = 0 ; x = 0 ; dx = 1 ; var = 0     
               else:
                  print "continue binning"
                  
            except IOError, (errno, strerror):
               print "I/O Error(%s) : (%s)" % (errno, strerror)      
            except:   
               print 'record skipped '+rec 
   f.close()

I have a script that round-trips correctly from a QDP file. I use it to read data from Swift XRT light curves produced by their online page.
Adapting to the Table reader/writer interface should not be too complicated, I'll look into it

The main difficulty is that a QDP file often contains more than one table, so it is not straightforward to adapt to the usual Table.read() interface.
Also, it has no standard header with column names. For example, Swift's header contains a command with column names, but adds a space in the middle of the Err (pos) column name, but not Err(neg), making it non-machine readable (with standard Astropy readers, at least):

READ TERR 1
READ SERR 2
! WT -- soft data
!MJD            Err (pos)       Err(neg)        Rate            Error
53000.123456 2.37847222222222e-05    -2.37847222222222e-05   0.726155        0.583890

So, what I did is something that interprets correctly the commands specifying the error columns (e.g. READ SERR 1, READ TERR 2) and gives some standard column names, ignoring the column names in the header. Then, returns a list of tables.

See here: https://gist.github.com/matteobachetti/5d102ce518bd183bf48b0aec2b4fa793

For reading, the table could be selected with index=<table_number> in the read() call.

I like Tom's proposal. When I read a QDP file, I do that anyway a table at a time, and this would fit better with the concept that one table is returned in a call. Having sometimes return a Table object and at other times a list means testing for that, and I do not see why the overhead. Rather than index=<table_number> use "table_number=<table_number>" which is more descriptive. The word index is overused. Paul

…

On Sat, Jan 16, 2021 at 11:11 PM Tom Aldcroft ***@***.***> wrote: For reading, the table could be selected with index=<table_number> in the read() call. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#8382 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABHLWCDGU64BGVRQ7SXGLC3S2IMILANCNFSM4GSRCYHA> .

-- * * * * * * * * http://www.mssl.ucl.ac.uk/~npmk/ * * * * N.P.M. Kuin ([email protected]) -- phone +44-(0)1483 (prefix) -204111 (work) -- mobile +44(0)7908715953 -- skype ID: npkuin Mullard Space Science Laboratory – University College London – Holmbury St Mary – Dorking – Surrey RH5 6NT– U.K.

Navigation Menu

Repository navigation

IO support for reading/writing QDP files #8382

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

IO support for reading/writing QDP files #8382

Description

Activity

taldcroft commented on Jan 26, 2019

PaulKuin commented on Jan 26, 2019

taldcroft commented on Jan 26, 2019

PaulKuin commented on Jan 26, 2019

PaulKuin commented on Jan 26, 2019

matteobachetti commented on Jan 15, 2021

matteobachetti commented on Jan 15, 2021

taldcroft commented on Jan 16, 2021

PaulKuin commented on Jan 17, 2021

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions