new release 1.1, simpler plot

GAIA_query.py

@@ -53,14 +53,19 @@ def main():
             print("Filters applied, {:.1f}% of data lost".format(
                 100. - (len(data) * 100.) / N_old))
 
-        print("Obtaining magnitudes/colors and their uncertainties")
-        mag, e_mag, col1, e_col1, col2, e_col2 = uncertMags(
-            data, col1_n, col2_n)
-
-        print("Write output file in input/ folder")
-        writeOut.main(
-            clust['name'], data, mag, e_mag, col1, e_col1, col2, e_col2,
-            col1_n, col2_n)
+        if read is False:
+            print("Obtaining magnitudes/colors and their uncertainties")
+            mag, e_mag, col1, e_col1, col2, e_col2 = uncertMags(
+                data, col1_n, col2_n)
+
+            print("Write output file in input/ folder")
+            writeOut.main(
+                clust['name'], data, mag, e_mag, col1, e_col1, col2, e_col2,
+                col1_n, col2_n)
+        else:
+            mag, e_mag, col1, e_col1, col2, e_col2 = data['Gmag'],\
+                data['e_Gmag'], data[col1_n], data['e_' + col1_n],\
+                data[col2_n], data['e_' + col2_n]
 
         print("Plotting")
         makePlot.main(
@@ -71,7 +76,9 @@ def main():
             mag, e_mag, col1, e_col1, col2, e_col2,
             data['Plx'], data['pmRA'], data['e_pmRA'], data['pmDE'],
             data['e_pmDE'], data['RV'], col1_n, col2_n,
-            babusiaux_filters)
+            babusiaux_filters, clust['box_s'])
+
+    print("\nEnd")
 
 
 def babusiaux_filt(data):

modules/_version.py

@@ -1 +1 @@
-__version__ = "v1.0"
+__version__ = "v1.1"

modules/makePlot.py

@@ -3,19 +3,18 @@
 from scipy.spatial import distance
 import matplotlib.pyplot as plt
 import matplotlib.cm as cm
-from matplotlib.colors import Normalize
 import matplotlib.offsetbox as offsetbox
 import matplotlib.gridspec as gridspec
 import warnings
 
 
 def main(
-    name, center, rad, e_mmax, e_c1max, e_c2max, plx_min, plx_max,
-        ra, dec, mag, e_mag, col1, e_col1, col2, e_col2,
-        plx, pmRA, e_pmRA, pmDE, e_pmDE, radv, col1_n, col2_n,
-        babusiaux_filters):
+    name, center, rad, e_mmax, e_c1max, e_c2max, plx_min, plx_max, ra, dec,
+    mag, e_mag, col1, e_col1, col2, e_col2, plx, pmRA, e_pmRA, pmDE, e_pmDE,
+        radv, col1_n, col2_n, babusiaux_filters, box_s, Nmax=25000):
     """
     """
+
     plt.style.use('seaborn-darkgrid')
     fig = plt.figure(figsize=(30, 25))
     gs = gridspec.GridSpec(10, 12)
@@ -26,28 +25,43 @@ def main(
     for i, col in enumerate([ra, mag, col1, col2, plx, pmRA, radv]):
         # Count valid data for each column
         col_sizes[0].append(col_names[i])
-        col_sizes[1].append(col[~col.mask].size)
+        try:
+            col_sizes[1].append(col[~col.mask].size)
+        except AttributeError:
+            col_sizes[1].append(col.size)
 
     ax1 = plt.subplot(gs[0:2, 0:2])
+    ax1.set_title("Field: [{} x {}]".format(box_s, box_s), fontsize=10)
     ax1.bar(col_sizes[0], col_sizes[1])
     fig.autofmt_xdate()
 
+    # Select at most Nmax random stars to plot
+    if ra.size > Nmax:
+        print("WARNING: too many stars ({}). Plotting {} random stars".format(
+            ra.size, Nmax))
+        ra, dec, mag, e_mag, col1, e_col1, col2, e_col2, plx, pmRA, e_pmRA,\
+            pmDE, e_pmDE, radv = random_choice(
+                ra, dec, mag, e_mag, col1, e_col1, col2, e_col2, plx, pmRA,
+                e_pmRA, pmDE, e_pmDE, radv, Nmax)
+
     ax2 = plt.subplot(gs[0:2, 2:4])
     ax2.set_title(r"$N_{{T}}={},\, rad={:.3f}\,[deg]$".format(
-        ra.size, rad), fontsize=8)
+        ra.size, rad), fontsize=10)
     plt.xlabel("RA [deg]")
     plt.ylabel("DEC [deg]")
     ra_min = max(min(ra), center[0] - 3. * rad)
     ra_max = min(max(ra), center[0] + 3. * rad)
+    ra_rng = ra_max - ra_min
     de_min = max(min(dec), center[1] - 3. * rad)
     de_max = min(max(dec), center[1] + 3. * rad)
-    msk = (ra > ra_min) & (ra < ra_max) & (dec > de_min) & (dec < de_max)
-    ax2.scatter(ra[msk], dec[msk], s=star_size(mag[msk]), c='k')
+    de_rng = de_max - de_min
+    rade_rng = min(ra_rng, de_rng) / 2.
+    ax2.scatter(ra, dec, s=star_size(mag), c='k')
     # Radius
     circle = plt.Circle(center, rad, color='red', lw=1.5, fill=False)
     fig.gca().add_artist(circle)
-    plt.xlim(ra_min, ra_max)
-    plt.ylim(de_min, de_max)
+    plt.xlim(center[0] - rade_rng, center[0] + rade_rng)
+    plt.ylim(center[1] - rade_rng, center[1] + rade_rng)
     ax2.invert_xaxis()
 
     ax4 = plt.subplot(gs[0:2, 4:6])
@@ -71,44 +85,38 @@ def main(
     ax5 = plt.subplot(gs[2:4, 0:2])
     ax5.set_title(
         "N(r<{}, e_mag<{}, e_c1<{}, e_c2<{})={}".format(
-            rad, e_mmax, e_c1max, e_c2max, msk.data.sum()), fontsize=8)
+            rad, e_mmax, e_c1max, e_c2max, msk.data.sum()), fontsize=10)
     plt.xlabel(col1_n)
     plt.ylabel(mag_n)
     ax5.scatter(col1[msk], mag[msk], s=4, lw=.1, edgecolor='w')
-    # no masked elements
-    msk4 = (~col1[msk].mask) & (~mag[msk].mask)
     x_max_cmd, x_min_cmd, y_min_cmd, y_max_cmd = diag_limits(
-        'mag', col1[msk][msk4], mag[msk][msk4])
+        'mag', col1[msk], mag[msk])
     plt.xlim(x_min_cmd, x_max_cmd)
     plt.ylim(y_min_cmd, y_max_cmd)
 
     ax6 = plt.subplot(gs[2:4, 2:4])
     plt.xlabel(col2_n)
     plt.ylabel(mag_n)
     ax6.scatter(col2[msk], mag[msk], s=4, lw=.1, edgecolor='w')
-    # no masked elements
-    msk4 = (~col2[msk].mask) & (~mag[msk].mask)
     x_max_cmd, x_min_cmd, y_min_cmd, y_max_cmd = diag_limits(
-        'mag', col2[msk][msk4], mag[msk][msk4])
+        'mag', col2[msk], mag[msk])
     plt.xlim(x_min_cmd, x_max_cmd)
     plt.ylim(y_min_cmd, y_max_cmd)
 
     ax9 = plt.subplot(gs[2:4, 4:6])
     plt.xlabel(col1_n)
     plt.ylabel(col2_n)
     ax9.scatter(col1[msk], col2[msk], s=4, lw=.1, edgecolor='w')
-    # no masked elements
-    msk4 = (~col1[msk].mask) & (~col2[msk].mask)
     x_max_cmd, x_min_cmd, y_min_cmd, y_max_cmd = diag_limits(
-        'col', col1[msk][msk4], col2[msk][msk4])
+        'col', col1[msk], col2[msk])
     plt.xlim(x_min_cmd, x_max_cmd)
     plt.ylim(y_min_cmd, y_max_cmd)
 
     if plx_min < 0.:
         print("Setting plx_min=0.")
         plx_min = 0.
     ax8 = plt.subplot(gs[4:6, 0:2])
-    ax8.set_title("{} < Plx [mas] < {}".format(plx_min, plx_max), fontsize=8)
+    ax8.set_title("{} < Plx [mas] < {}".format(plx_min, plx_max), fontsize=10)
     plt.xlabel("Plx [mas]")
     msk2 = (d_col < rad) & (e_mag < e_mmax) &\
         (e_col1 < e_c1max) & (e_col2 < e_c2max) &\
@@ -126,38 +134,26 @@ def main(
     ax8.add_artist(ob)
     plt.xlim(0., 3.)
 
-    ax7 = plt.subplot(gs[4:6, 2:4])
-    ax7.set_title("N(r<, e_G<, e_XP<, <Plx<)={}".format(
-        msk2.data.sum()), fontsize=8)
-    plt.xlabel("RA [deg]")
-    plt.ylabel("DEC [deg]")
-    cmap = cm.viridis_r
-    norm = Normalize(vmin=0., vmax=p_max_mas)
-    ax7.scatter(ra[msk2], dec[msk2], s=4, c=cmap(norm(plx[msk2])))
-    ax7.invert_xaxis()
-    im = plt.scatter(ra[msk2], dec[msk2], s=0, c=plx[msk2], cmap=cmap)
-    cbar_ax = fig.add_axes([0.313, 0.53, 0.005, 0.05])
-    cbar = fig.colorbar(im, cax=cbar_ax)
-    cbar.ax.tick_params(labelsize=5)
-    plt.clim(0., 2. * p_max_mas)
-
-    ax3 = plt.subplot(gs[4:6, 4:6])
+    ax3 = plt.subplot(gs[4:6, 2:4])
     plt.xlabel("pmRA [mas/yr]")
     plt.ylabel("pmDEC [mas/yr]")
-    msk3 = (pmRA > -30) & (pmDE > -30) & (pmRA < 30) &\
-        (pmDE < 30) & (d_col < rad)
-    pmRA_f, pmDE_f, epmRA_f, epmDE_f = pmRA[msk3], pmDE[msk3],\
-        e_pmRA[msk3], e_pmDE[msk3]
-    ax3.set_title("N(r<rad, |pmX|<30)={}".format(pmRA_f.size), fontsize=8)
-    cmap = cm.viridis
-    norm = Normalize(vmin=d_col[msk3].min(), vmax=rad)
-    ax3.errorbar(
-        pmRA_f, pmDE_f, yerr=epmDE_f, xerr=epmRA_f, fmt='none', elinewidth=.35,
-        ecolor=cmap(norm(d_col[msk3])))
-    im = plt.scatter(pmRA_f, pmDE_f, s=0, c=d_col[msk3], cmap=cmap)
-    cbar_ax = fig.add_axes([0.48, 0.53, 0.005, 0.05])
-    cbar = fig.colorbar(im, cax=cbar_ax)
-    cbar.ax.tick_params(labelsize=5)
+    msk_in = d_col <= rad
+    rain_mean, rain_std = pmRA[msk_in].mean(), pmRA[msk_in].std()
+    dein_mean, dein_std = pmDE[msk_in].mean(), pmDE[msk_in].std()
+    ra_min, ra_max = rain_mean - 5. * rain_std, rain_mean + 5. * rain_std
+    de_min, de_max = dein_mean - 5. * dein_std, dein_mean + 5. * dein_std
+    ax3.set_title(
+        "N(r<rad, |pmX|<30)={}".format(pmRA[msk_in].size), fontsize=10)
+    ax3.plot(
+        pmRA[~msk_in], pmDE[~msk_in], 'k+', alpha=.25, zorder=1)
+    ax3.scatter(
+        pmRA[msk_in], pmDE[msk_in], c=d_col[msk_in], s=10, alpha=.25,
+        cmap=cm.viridis, zorder=5)
+    plt.xlim(ra_min, ra_max)
+    plt.ylim(de_min, de_max)
+    ax3.invert_xaxis()
+
+    # ax7 = plt.subplot(gs[4:6, 4:6])
 
     # Ignore warning issued by colorbar.
     with warnings.catch_warnings():
@@ -189,12 +185,9 @@ def diag_limits(yaxis, phot_x, phot_y):
     '''
     Define plot limits for *all* photometric diagrams.
     '''
-    # TODO deprecated
-    # min_x, max_x, min_y, max_y = kde_limits(phot_x, phot_y)
-
-    x_median, x_std = np.median(phot_x), 1.5 * np.std(phot_x)
+    x_median, x_std = np.ma.median(phot_x), 1.5 * np.ma.std(phot_x)
     min_x, max_x = x_median - x_std, x_median + x_std
-    y_median, y_std = np.median(phot_y), np.std(phot_y)
+    y_median, y_std = np.ma.median(phot_y), np.ma.std(phot_y)
     min_y, max_y = y_median - y_std, y_median + y_std
 
     # Define diagram limits.
@@ -208,3 +201,15 @@ def diag_limits(yaxis, phot_x, phot_y):
         y_max_cmd = min_y - 1.
 
     return x_max_cmd, x_min_cmd, y_min_cmd, y_max_cmd
+
+
+def random_choice(
+    ra, dec, mag, e_mag, col1, e_col1, col2, e_col2, plx, pmRA, e_pmRA, pmDE,
+        e_pmDE, radv, Nmax):
+    """
+    """
+    idxs = np.random.choice(ra.size, Nmax)
+
+    return ra[idxs], dec[idxs], mag[idxs], e_mag[idxs], col1[idxs],\
+        e_col1[idxs], col2[idxs], e_col2[idxs], plx[idxs], pmRA[idxs],\
+        e_pmRA[idxs], pmDE[idxs], e_pmDE[idxs], radv[idxs]