I have a catalogue of 56000 objects and I want to replace some columns with new values for instance replace NaN with -99 and 0 for specific columns (as you can see in the following code). Then I need to compute a value for some columns with a function which is called  get_ebv_from_mapand subtract it from some other columns and so on.Using get_ebv_from_map function to estimate the values in the catalogue is taking too much time and I think what I did like using this loop is not very efficient.The function is given by:from astropy import units as ufrom astropy.coordinates import SkyCoordfrom scipy.ndimage.interpolation import map_coordinatesimport os.pathfrom collections import Iterableimport pyfits as fitsdef get_ebv_from_map(coordinates, interpolate=True, order=1):    mapdir='/vol/Dust_Map/maps/'    fname = os.path.join(mapdir, 'SFD_dust_4096_{0}.fits')    # Parse input     ra, dec = coordinates    coordinates = SkyCoord(ra=ra, dec=dec, unit=(u.degree, u.degree))    # Convert to galactic coordinates.    coordinates = coordinates.galactic    l = coordinates.l.radian    b = coordinates.b.radian    # Check if l, b are scalar    return_scalar = False    if not isinstance(l, Iterable):        return_scalar = True        l, b = np.array([l]), np.array([b])    # Initialize return array    ebv = np.empty_like(l)    # Treat north (b>0) separately from south (b<0).    for n, idx, ext in [(1, b >= 0, 'ngp'), (-1, b < 0, 'sgp')]:        if not np.any(idx): continue        hdulist = fits.open(fname.format(ext))        mapd = hdulist[0].data        # Project from galactic longitude/latitude to lambert pixels.        # (See SFD98).        npix = mapd.shape[0]                x = (npix / 2 * np.cos(l[idx]) * np.sqrt(1. - n*np.sin(b[idx])) +             npix / 2 - 0.5)        y = (-npix / 2 * n * np.sin(l[idx]) * np.sqrt(1. - n*np.sin(b[idx])) +             npix / 2 - 0.5)        # Get map values at these pixel coordinates.        if interpolate:            ebv[idx] = map_coordinates(mapd, [y, x], order=order)        else:            x=np.round(x).astype(np.int)            y=np.round(y).astype(np.int)            ebv[idx] = mapd[y, x]        hdulist.close()    if return_scalar:        return ebv[0]    return ebvThe loop is given in below:ABtoVEGA=[+0.77,-0.13,-0.02,+0.19,+0.49,-0.19,-0.18,-0.06,-0.06,+0.04,+0.10,+0.22,+0.27,+0.36,+0.45,+0.56,+0.50]ZP=np.loadtxt(zero-points.offsets.asc)Gal=np.loadtxt(photometry_galaxies_photometry.cat, dtype=None)ra=Gal[:,1]dec=Gal[:,2]Av=[3.836 ,3.070,2.542,2.058,1.313, 3.466 , 3.082 , 2.882 , 2.652 , 2.367 , 2.226 , 2.066 , 1.872 , 1.652 , 1.423 , 1.298 , 1.156]count=0for i in range(Gal.shape[0]):    j=0    m=39    while (j<len(ABtoVEGA)):         if ((isnan(float(Gal[i,j+5]))) or (isnan(float(Gal[i,j+22])))):            Gal[i,j+5]='-99.0'            Gal[i,j+22]='0.0'            count+=1         else:            value=get_ebv_from_map((ra[i], dec[i]), interpolate=True, order=1)            Gal[i,j+5]=str(float(Gal[i,j+5])-ABtoVEGA[j]-ZP[j]-value*Av[j])         j+=1    while (m<45):          if (math.isnan(float( Gal[i,m]))):              Gal[i,m] = '0.0'          m+=1Do you have any advice on how I could make it optimal and fast?
Catalogue manipulation
python;vectorization;multiprocessing;scipy;astropy