In [2]:
import numpy as np
from astropy.table import Table
from wquantiles import median as wmedian
import os
In [3]:
from mltier1 import describe, Field
In [4]:
%pylab inline
This is the object representing the region considered. It can be used to filter the catalogues and to compute the area.
In [5]:
field = Field(170.0, 190.0, 45.5, 56.5)
In [6]:
field_full = Field(160.0, 232.0, 42.0, 62.0)
In [8]:
base_panstarrs = Table.read(os.path.join("data", "hetdex_full_jsabater.fits"))
In [9]:
len(base_panstarrs)
Out[9]:
The original PanSTARRS catalogue had ~4.4M sources but the catalogue that covers all the area has ~27M sources. The columns are the following:
In [10]:
np.array(base_panstarrs.colnames)
Out[10]:
In [11]:
describe(base_panstarrs['decMean'])
describe(base_panstarrs['raMean'])
describe(base_panstarrs['decStack'])
describe(base_panstarrs['raStack'])
In [12]:
describe(base_panstarrs['iFApFlux'])
describe(base_panstarrs['iFApFluxErr'])
In [13]:
hexbin(base_panstarrs["raMean"], base_panstarrs["decMean"]);
In [14]:
plot(base_panstarrs["raMean"][np.isnan(base_panstarrs['iFApFlux'])],
base_panstarrs["decMean"][np.isnan(base_panstarrs['iFApFlux'])],
ls="", marker=",");
The next figure is the distribution of sources with wrong values for the forced photometry in i-band but which were measured.
In [15]:
plot(base_panstarrs["raMean"][base_panstarrs['iFApFlux'] == -999.],
base_panstarrs["decMean"][base_panstarrs['iFApFlux'] == -999.],
ls="", marker=",");
In [16]:
offset_dec = base_panstarrs['decMean']-base_panstarrs['decStack']
offset_ra = base_panstarrs['raMean']-base_panstarrs['raStack']
In [17]:
describe(offset_ra, decimals=7)
describe(offset_dec, decimals=7)
In [18]:
hexbin(offset_ra, offset_dec, bins='log')
Out[18]:
In [19]:
describe(base_panstarrs['gPSFMag'])
describe(base_panstarrs['gPSFMagErr'])
In [20]:
describe(base_panstarrs['rPSFMag'])
describe(base_panstarrs['rPSFMagErr'])
In [21]:
describe(base_panstarrs['iPSFMag'])
describe(base_panstarrs['iPSFMagErr'])
In [22]:
describe(base_panstarrs['zPSFMag'])
describe(base_panstarrs['zPSFMagErr'])
In [23]:
describe(base_panstarrs['yPSFMag'])
describe(base_panstarrs['yPSFMagErr'])
In [24]:
len(base_panstarrs[
((base_panstarrs["gPSFMag"] == - 999.) &
(base_panstarrs["rPSFMag"] == - 999.) &
(base_panstarrs["iPSFMag"] == - 999.) &
(base_panstarrs["zPSFMag"] == - 999.) &
(base_panstarrs["yPSFMag"] == - 999.))])
Out[24]:
In [25]:
new_panstarrs = field_full.filter_catalogue(
base_panstarrs[
((base_panstarrs["gPSFMag"] != - 999.) |
(base_panstarrs["rPSFMag"] != - 999.) |
(base_panstarrs["iPSFMag"] != - 999.) |
(base_panstarrs["zPSFMag"] != - 999.) |
(base_panstarrs["yPSFMag"] != - 999.))],
colnames=("raMean", "decMean"))
In [26]:
len(new_panstarrs)
Out[26]:
Free memory
In [27]:
del base_panstarrs
In [29]:
panstarrs = new_panstarrs.to_pandas()
In the auxiliary catalogue we will use NaN to define the missing values for the magnitudes.
In [30]:
all_bands = ["g", "r", "i", "z", "y"]
In [31]:
for band in all_bands:
panstarrs.loc[(panstarrs["%sPSFMag"%band] == -999.), "%sPSFMag"%band] = np.nan
panstarrs.loc[(panstarrs["%sPSFMagErr"%band] == -999.), "%sPSFMagErr"%band] = np.nan
In [32]:
len(panstarrs[np.isnan(panstarrs["iPSFMag"])])
Out[32]:
We get the statistical properties for each of the other bands with data
In [33]:
bands = ["g", "r", "z", "y"]
In [34]:
def get_diff(band):
assert band in bands
column = "{}PSFMag".format(band)
mask = ((panstarrs[column] != -999.) &
(panstarrs["iPSFMag"] != -999.) &
(~np.isnan(panstarrs[column])) &
(~np.isnan(panstarrs["iPSFMag"])))
return panstarrs[mask][column] - panstarrs[mask]["iPSFMag"]
In [35]:
band_properties = {}
for band in bands:
diff = get_diff(band)
band_properties.update({band: {"mean": diff.mean(),
"median": np.median(diff),
"std": diff.std(),
"n": len(diff)}
})
if True:
print("{} {n} {mean:7.4f} {std:7.4f} {median:7.4f}".format(band, **band_properties[band]))
Now we compute the interpolated value of "i" for the missing data. The computed values, including the measured ones, are entered in a column named "i".
In [36]:
aux_weights = np.array([1,2,2,1])
mags_i_diff = np.array([band_properties[a]["mean"] for a in bands])
In [37]:
def interpolate_i(x):
if np.isnan(x["iPSFMag"]):
mags_raw = x[["gPSFMag",
"rPSFMag",
"zPSFMag",
"yPSFMag"]].values
mags_i = mags_raw - mags_i_diff
try:
return wmedian(mags_i[~np.isnan(mags_raw)],
aux_weights[~np.isnan(mags_raw)])
except ValueError:
return np.nan
else:
return x["iPSFMag"]
In [38]:
panstarrs["i"] = panstarrs.apply(interpolate_i, axis=1)
The error in the measurement is also computed and entered in a new column.
In [39]:
aux_err = np.array([band_properties[a]["std"]/np.sqrt(band_properties[a]["n"]) for a in bands])
def interpolate_error_i(x):
if np.isnan(x["iPSFMagErr"]):
err_mags_raw = x[["gPSFMagErr",
"rPSFMagErr",
"zPSFMagErr",
"yPSFMagErr"]].values
errs1 = err_mags_raw[~np.isnan(err_mags_raw)]
errs2 = aux_err[~np.isnan(err_mags_raw)]
try:
return np.sqrt(np.sum(errs1**2) + np.sum(errs2**2))/np.sqrt(0.5*(len(errs1)+len(errs2)))
except ValueError:
return np.nan # There are some sources here with no errors in any band FIXME
else:
return x["iPSFMagErr"]
In [40]:
panstarrs["iErr"] = panstarrs.apply(interpolate_error_i, axis=1)
In [41]:
panstarrs["i"].describe()
Out[41]:
In [42]:
panstarrs["iErr"].describe()
Out[42]:
In [43]:
new_panstarrs["i"] = panstarrs["i"].values
new_panstarrs["iErr"] = panstarrs["iErr"].values
In [44]:
columns_save = ['objID', 'raMean', 'decMean', 'raMeanErr', 'decMeanErr', 'i', 'iErr']
In [45]:
new_panstarrs[columns_save].write('panstarrs_u1.fits', format="fits")
Free memory
In [46]:
del panstarrs
del new_panstarrs
In [47]:
base_wise = Table.read(os.path.join("data", "allwise_HETDEX_full.fits"))
In [49]:
len(base_wise)
Out[49]:
In [50]:
np.array(base_wise.colnames)
Out[50]:
In [51]:
describe(base_wise['ra'])
describe(base_wise['dec'])
In [52]:
hexbin(base_wise["ra"], base_wise["dec"]);
In [53]:
new_wise = field_full.filter_catalogue(base_wise)
In [54]:
new_wise["W1mag"] = new_wise['w1mpro'] + 2.683
In [55]:
new_wise.rename_column("designation", 'AllWISE')
new_wise.rename_column("ra", 'raWise')
new_wise.rename_column("dec", 'decWise')
new_wise.rename_column("sigra", 'raWiseErr')
new_wise.rename_column("sigdec", 'decWiseErr')
new_wise.rename_column('w1sigmpro', "W1magErr")
In [56]:
columns_save = ['AllWISE', 'raWise', 'decWise', 'raWiseErr', 'decWiseErr', 'W1mag', 'W1magErr']
In [ ]:
new_wise[columns_save].write('wise_u1.fits', format="fits")
In [ ]: