diff --git a/stellarphot/notebooks/photometry/04-transform-pared-back.ipynb b/stellarphot/notebooks/photometry/04-transform-pared-back.ipynb old mode 100755 new mode 100644 index fef0dbd7..a847b467 --- a/stellarphot/notebooks/photometry/04-transform-pared-back.ipynb +++ b/stellarphot/notebooks/photometry/04-transform-pared-back.ipynb @@ -19,11 +19,13 @@ "import matplotlib.pyplot as plt\n", "\n", "from stellarphot.utils.magnitude_transforms import (\n", - " f,\n", + " calibrated_from_instrumental,\n", " opts_to_str,\n", " calc_residual,\n", " transform_to_catalog,\n", - ")" + ")\n", + "\n", + "from stellarphot import PhotometryData, vsx_vizier" ] }, { @@ -74,29 +76,49 @@ "c_min = -0.5\n", "d_min = -1e-6\n", "\n", - "our_filters = [\"SI\"]\n", + "our_filters = [\"SR\"]\n", + "\n", + "aavso_band_names = dict(B=\"B\", V=\"V\", gp=\"SG\", rp=\"SR\", ip=\"SI\", SI=\"SI\", SR=\"SR\")\n", "\n", - "aavso_band_names = dict(B=\"B\", V=\"V\", gp=\"SG\", rp=\"SR\", ip=\"SI\", SI=\"SI\")\n", + "# cat_color_colums = dict(\n", + "# B=(\"Bmag\", \"Vmag\"),\n", + "# V=(\"Bmag\", \"Vmag\"),\n", + "# gp=(\"g_mag\", \"r_mag\"),\n", + "# rp=(\"r_mag\", \"i_mag\"),\n", + "# ip=(\"r_mag\", \"i_mag\"),\n", + "# SI=(\"r_mag\", \"i_mag\"),\n", + "# )\n", "\n", "cat_color_colums = dict(\n", - " B=(\"Bmag\", \"Vmag\"),\n", - " V=(\"Bmag\", \"Vmag\"),\n", - " gp=(\"g_mag\", \"r_mag\"),\n", - " rp=(\"r_mag\", \"i_mag\"),\n", - " ip=(\"r_mag\", \"i_mag\"),\n", - " SI=(\"r_mag\", \"i_mag\"),\n", + " B=(\"mag_B\", \"mag_V\"),\n", + " V=(\"mag_B\", \"mag_V\"),\n", + " gp=(\"mag_SG\", \"mag_SR\"),\n", + " rp=(\"mag_SR\", \"mag_SI\"),\n", + " SR=(\"SR\", \"SI\"),\n", + " ip=(\"mag_SR\", \"mag_SI\"),\n", + " SI=(\"mag_SR\", \"mag_SI\"),\n", ")\n", "\n", + "# cat_filter = dict(\n", + "# B=\"Bmag\",\n", + "# V=\"Vmag\",\n", + "# gp=\"g_mag\",\n", + "# rp=\"r_mag\",\n", + "# ip=\"i_mag\",\n", + "# SI=\"i_mag\",\n", + "# )\n", + "\n", "cat_filter = dict(\n", - " B=\"Bmag\",\n", - " V=\"Vmag\",\n", - " gp=\"g_mag\",\n", - " rp=\"r_mag\",\n", - " ip=\"i_mag\",\n", - " SI=\"i_mag\",\n", + " B=\"B\",\n", + " V=\"mag_V\",\n", + " gp=\"mag_SG\",\n", + " rp=\"mag_SR\",\n", + " SR=\"SR\",\n", + " ip=\"mag_SI\",\n", + " SI=\"mag_SI\",\n", ")\n", "\n", - "input_photometry_file = \"TIC-402828941-2021-09-28.ecsv\"\n", + "input_photometry_file = \"sp2-kelt-1b-photometry-2021-09-14.ecsv\"\n", "output_photometry_file = \"some_name.ecsv\"" ] }, @@ -108,7 +130,11 @@ }, "outputs": [], "source": [ - "all_mags = Table.read(input_photometry_file)\n", + "all_mags = PhotometryData.read(input_photometry_file)\n", + "\n", + "# BAD BAD BAD BAD BAD BAD BAD BAD BAD BAD\n", + "all_mags[\"passband\"] = \"SR\"\n", + "\n", "\n", "# Ensure we have the right table ordering later\n", "all_mags.sort([\"passband\", \"bjd\"])" @@ -137,6 +163,15 @@ "assert all(k[0] in cat_filter.keys() for k in filter_groups.groups.keys)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "filter_groups.groups.keys" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -148,6 +183,7 @@ "cell_type": "code", "execution_count": null, "metadata": { + "scrolled": true, "tags": [] }, "outputs": [], @@ -156,7 +192,7 @@ "\n", "for k, group in zip(our_filters, filter_groups.groups):\n", " print(f\"Transforming band {k}\")\n", - " by_bjd = group.group_by(\"bjd\")\n", + " by_bjd = group.group_by(\"file\")\n", "\n", " transform_to_catalog(\n", " by_bjd,\n", @@ -190,7 +226,32 @@ "metadata": {}, "outputs": [], "source": [ - "plt.plot(by_bjd[\"mag_cat\"] - by_bjd[\"mag_inst_cal\"], \"o\")\n", + "plt.plot(by_bjd[\"mag_cat\"], by_bjd[\"mag_cat\"] - by_bjd[\"mag_inst_cal\"], \"o\")\n", + "plt.ylim(-0.05, 0.05)\n", + "plt.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "kelt1 = by_bjd.lightcurve_for(\"kelt-1\", flux_column=\"mag_inst_cal\") # by_bjd[\"star_id\"] == 1\n", + "plt.plot(kelt1.time.jd, kelt1.flux, \".\")\n", + "plt.ylim(reversed(plt.ylim()))\n", + "plt.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.plot(kelt1.time.value, kelt1[\"z\"], \".\")\n", "plt.grid()" ] }, @@ -203,6 +264,24 @@ "output_table.write(output_photometry_file)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "by_bjd.colnames" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "by_bjd['id']" + ] + }, { "cell_type": "code", "execution_count": null, @@ -227,7 +306,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.11.9" } }, "nbformat": 4, diff --git a/stellarphot/utils/magnitude_transforms.py b/stellarphot/utils/magnitude_transforms.py index 53d864f9..f4c1c719 100644 --- a/stellarphot/utils/magnitude_transforms.py +++ b/stellarphot/utils/magnitude_transforms.py @@ -424,7 +424,7 @@ def transform_magnitudes( transform_catalog_coords = SkyCoord( transform_catalog["RAJ2000"], transform_catalog["DEJ2000"], unit="deg" ) - input_coords = SkyCoord(input_mags["RA"], input_mags["Dec"]) + input_coords = SkyCoord(input_mags["ra"], input_mags["dec"]) transform_catalog_index, d2d, _ = match_coordinates_sky( input_coords, transform_catalog_coords @@ -593,22 +593,51 @@ def transform_to_catalog( cat = None cat_coords = None + def fake_it(one_image): + # Accumulate the parameters + popt = [np.nan] * len(all_params) + for param, value in zip(all_params, popt, strict=True): + param.extend([value] * len(one_image)) + nana = [np.nan] * len(one_image) + cal_mags.extend(nana) + cat_mags.extend(nana) + cat_colors.extend(nana) + resids.extend(nana) + for file, one_image_inp in zip( observed_mags_grouped.groups.keys, observed_mags_grouped.groups, strict=True ): one_image = one_image_inp[one_image_inp["passband"] == obs_filter] - our_coords = SkyCoord(one_image["RA"], one_image["Dec"], unit="degree") + our_coords = SkyCoord(one_image["ra"], one_image["dec"], unit="degree") if cat is None or cat_coords is None: cat = apass_dr9( our_coords[0], radius=1 * u.degree, clip_by_frame=False, padding=0 ) - cat_coords = SkyCoord(cat["ra"], cat["dec"]) - cat["color"] = cat[cat_color[0]] - cat[cat_color[1]] + cat = cat.passband_columns(passbands=["B", "V", "SR", "SG", "SI"]) + cat["mag_RC"] = filter_transform( + cat, + output_filter="R", + g="mag_SG", + r="mag_SR", + i="mag_SI", + transform="jester", + ) + cat["mag_IC"] = filter_transform( + cat, + output_filter="I", + g="mag_SG", + r="mag_SR", + i="mag_SI", + transform="jester", + ) + + cat_coords = SkyCoord(cat["ra"], cat["dec"], unit="degree") + cat["color"] = cat[f"mag_{cat_color[0]}"] - cat[f"mag_{cat_color[1]}"] cat_idx, d2d, _ = our_coords.match_to_catalog_sky(cat_coords) mag_inst = one_image[obs_mag_col] - cat_mag = cat[cat_filter][cat_idx] + cat_mag = cat[f"mag_{cat_filter}"][cat_idx] color = cat["color"][cat_idx] # Impose some constraints on what is included in the fit @@ -619,8 +648,22 @@ def transform_to_catalog( & ~np.isnan(one_image[obs_mag_col]) ) + if not (any(good_dat) or any(good_cat)): + print(f"No good data in {file[0]}") + fake_it(one_image) + continue + mag_diff = cat_mag - mag_inst + if ( + (not (any(good_dat) or any(good_cat))) + or all(np.isnan(mag_diff)) + or all(mag_diff.mask) + ): + print(f"No good data in {file[0]}") + fake_it(one_image) + continue + good_data = good_dat & ( np.abs(mag_diff - np.nanmedian(mag_diff[good_dat & ~mag_diff.mask])) < 1 ) @@ -631,6 +674,11 @@ def transform_to_catalog( good = good_cat & good_data + if not any(good): + print(f"No good data in {file[0]}") + fake_it(one_image) + continue + # Prep for fitting init_guess = (a_cen, 0, 0, 0, zero_point_mid) X = (mag_inst[good], color[good]) diff --git a/stellarphot/utils/tests/test_magnitude_transforms.py b/stellarphot/utils/tests/test_magnitude_transforms.py index b5eb95a4..65f27254 100644 --- a/stellarphot/utils/tests/test_magnitude_transforms.py +++ b/stellarphot/utils/tests/test_magnitude_transforms.py @@ -73,8 +73,8 @@ def generate_tables(n_stars, mag_model): ra, dec = generate_star_coordinates(n_stars) # Instrumental magnitudes - ra_col = Column(name="RA", data=ra) - dec_col = Column(name="Dec", data=dec) + ra_col = Column(name="ra", data=ra) + dec_col = Column(name="dec", data=dec) instrumental = Table([instr_mags, ra_col, dec_col])