Astrophysics Research Python Programming

Loading FITS and MATCH

The first we have to be able to read the ROTSE-I and ROTSE-III data files. There were three types of files.

1. .dat: MATCH structure
2. .datc: compressed MATCH structure
3. .fit: FITS structure

Luckily, all these structures are readable natively by python packages. The following code insert shows how to read MATCH and FITS structured files.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

from scipy.io import readsav import pyfits def read_fits_file(file, fits_index=1): try: hdus = pyfits.open(file, memmap=True) hdus_ext=hdus[fits_index] data=hdus_ext.data except Exception as e: raise Exception("cannot read fits data from file: %s" % (file,)) from e return data def read_match_file(file, *args, **kwargs): try: data=readsav(file)['match'] except Exception as e: raise Exception("cannot read match data from file: %s" % (file,)) from e return data

The above code uses two packages scipy.io and pyfits, which provides the actual loading of the structures. Note that readsav reads both .dat and .datc files.

The type of the resulting data is of type numpy.recarry. Elements of recarray can be accessed via name either by the '.' or '[]' notation.

It is easy to have a general purpose reader by having a routing function as follows:

1 2 3 4 5 6 7 8 9 10

def read_data_file(file, fits_index=1, tmpdir='/tmp'): if not os.path.isfile(file): raise Exception("file not found: %s" % (file,)) file_ext=file.rpartition('.')[2] if file_ext=='fit': data=read_fits_file(file, fits_index) else: # assume MATCH data=read_match_file(file) return data

Obviously, this is good only if the next processing steps takes into account differences between the two structures. For example, MATCH structure have a STAT field that is not part of FITS structute.

1 2	try: stat=data.field('STAT')[0] except: stat=None # assume not a MATCH file

Where statement

IDL programs take advantage of IDL's where_ statement. It will look something like follows:

1 2	select_indexes = where(data.flags[*,ptr] gt -1 and $ check_flags(ref_vflag, data.flags[*,ptr], type='STYPE') eq 0)

The above statement is pretty loaded in functionality. In essence, the expected result are indexes into the vector data.values[,k] that adheres to two conditions: the need to be greater than -1; and the result for the function *check_value to be zero.

1 2 3 4 5 6 7 8 9

import numpy as np # assum data is a MATCH or FITS data structure # and that check_flag function is defined for ref_flag and flags=data.field('FLAGS')[0] check_flags_v=np.vectorize(lambda flag: check_flag(ref_flag, flag) ==0) cond=np.logical_and.reduce( (flags[ptr,:] > -1, check_flags_v(flags[ptr,:]), ), ) select_indexes=np.where(cond); select_indexes=select_indexes[0]

To match IDL with Python, we are using four distinct numpy functions.

1. vectorize transforms a vector by applying a function on its elements.
2. logical_and and reduce: applys logical test to vector element transforming it to boolean vector.
3. where: return indices of those element set to True.

An important note is that IDL matrix indices are in opposite order to that of numpy.

From here on we will assume np stands for numpy imported as np.

Just to clarify, numpy performs operation on arrays. For example, assuming merr and msys are arrays of the same size, the following will produce a new array which each element is the square root of the squared sum of related elements. Obviously, a more sophisticated computation can be deployed.

1	np.sqrt(merr**2.0 + msys**2.0)

Hierarchical record array initialization

IDL's create_struct_ is being used to create records with fields accessible by name. It's parallel in numpy realm is recarray.

The following code set shows how to create recarray structure with default values. Starting with the definition of the records' fields, field_map. it is built as a list of tuple, each describe a name a field, a type, and a default value.

1 2 3 4 5 6 7

field_map=[ ('state', int, -1), ('distance', np.float32, -1.0), ('posangle', np.float32, -1.0), ('error', np.float32, -1.0), ('phot', np.float32, -1.0), ('photerror', np.float32, -1.0),]

Next, base on the above mapping, create_structured_vector would create the records' data type and a vector[size] with elements of the desired data type.

1 2 3 4 5 6 7 8 9 10 11

def create_structured_vector(size, field_map, copy=False): dts=list() for name, type_, _ in field_map: dts.append( (name, type_,) ) dtype = np.dtype(dts) values=[tuple([value if not copy else np.copy(value) for _, _, value in field_map]) for _ in range(size)] array=np.rec.array(values, dtype=dtype) return array, dtype

The copy option of create_structured_vector, if set, tells it to copy the default values. This is useful in case the default value is a structure by itself.

Note that the function returns a tuple of the generated array and the created type. This is useful in case further association of this data type is required.

This method could be extended to handle any array shape, not just vectors.

Plotting

IDL's plot capability can be achieved using matplotlib and pandas's dataframe.plot. These tools are rich with features and easy to use.

The original IDL code we were porting was built creating postscript documents. Using matplotlib we switched to PDF.

Here is an example how to print few drawing per page. We start with importing and setting matpolotlib for PDF plotting.

1 2 3 4

import matplotlib matplotlib.use('PDF') import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages

Next, we define subplot on a page organized by 3 rows of 2 figure each. Assuming that data is a list of tuples (title, values), we

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

pdf=PdfPages('pdffile.pdf') for i in range(len(data)): if i%6 == 0: # first in a page fig, axarr = plt.subplots(3, 2, figsize=(11, 8.5)) # data = [(title, x-values, y-values), ...] k=i%6; fig=axarr[int(k/2), k-int(k/2)*2] title, x_values, y_values = data[i] fig.set_title(title, fontsize=11) fig.plot(x_values, y_values, '+') if i%6 == 5: # last in a page fig.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0) pdf.savefig() plt.close() pdf.close()

Moment

moment is an IDL function that provides in a single call 4 statistical calculations. moment returns a four-element vector containing the mean, variance, skewness, and kurtosis of the input vector.

Similar functionality can be accomplish using numpy and scipy.

1
2
3
4
5
6
7
8
9

from scipy import stats

def moment(value_vec)
    mean=np.mean(value_vec)
    sdev=np.sqrt(np.var(value_vec))
    skew=stats.skew(value_vec)
    vkurt=stats.kurtosis(value_vec)

    return mean, sdev, skew, vkurt

scipy's stats_ has many other shortcut computations that are worth a while to look at. In reality, it is rare that one needs all four computational elements. Therefore, it is better to engage with the specific functions as needed instead of using moment().