Efficiently get midpoint of image in Python with astropy

My data is organized in a table in fits format (example.fits). I open this table ('mappars') using the Python module atpy. The table has columns x,y which are image coordinates and a data column z. The x,y are evenly spaced (events on a CCD detector) but there are gaps in between.

I am getting the value in the middle of the image like this:

import atpy
import numpy as np

mappars = atpy.Table('example.fits')
#get midpoint value
midx = np.int((np.max(mappars['x'])+np.min(mappars['x']))/2)
midy = np.int((np.max(mappars['y'])+np.min(mappars['y']))/2)
midist = mappars.where((mappars.x == midx) & (mappars.y == midy)['z']

Is there a more efficient way of doing this (without the .where function which is part of atpy)? Also is there a function corresponding to .where in astropy.table since I would like to move from atpy to astropy?

Solution

You should be able to do this very similarly using astropy.io and numpy. I don't have an event file handy, but something like:

from astropy.io import fits import numpy as np

with fits.open("example.fits") as hdulist:
    events = hdulist[1].data   # since hdu 0 should be an ImageHDU type, the events are probably in the first extension
    midx = np.int((np.max(events['x'])+np.min(events['x']))/2)
    midy = np.int((np.max(events['y'])+np.min(events['y']))/2)
    midist = events['z'][(events['x'] == midx) & (events['y'] == midy)]

In this case, I'm directly creating a boolean index array to events['z'] using the boolean (numpy) expression events['x'] == midx) & (events['y'] == midy).

Note that this is untested; if the above fails for you, I can probably try and pull an event file somewhere from the internet to test it, unless of course your event file is peculiar.

By the way, if x and y are, say, integer pixel coordinates, you can make things slightly simpler by using pure integer division and skip the cast to int:

midx = (np.max(events['x'])+np.min(events['x']))//2
midy = (np.max(events['y'])+np.min(events['y']))//2

This works in Python 3 directly, and in Python 2 if you have from __future__ import division (you probably should do that anyway).