Merge pull request #115 from astrophpeter/djones-prospector

Djones prospector

.gitignore

@@ -153,6 +153,7 @@ data/cutout_cdn/2010ai/*
 !data/cutout_cdn/2010h/*
 !data/cutout_cdn/2010ag/*
 !data/cutout_cdn/2010ai/*
+!data/sed_output/2010h/*
 
 static/
 

app/app/settings.py

@@ -134,6 +134,7 @@
 # os.path.join(os.path.dirname(BASE_DIR), '../cutout_cdn')
 
 MEDIA_ROOT = os.path.join(os.path.dirname(BASE_DIR), "../cutout_cdn")
+SED_OUTPUT_ROOT = os.path.join(os.path.dirname(BASE_DIR), "../sed_output")
 GHOST_OUTPUT_ROOT = os.path.join(os.path.dirname(BASE_DIR), "../ghost_output")
 TNS_STAGING_ROOT = os.path.join(os.path.dirname(BASE_DIR), "../tns_staging")
 TRANSMISSION_CURVES_ROOT = os.path.join(os.path.dirname(BASE_DIR), "../transmission")

app/entrypoints/clean_data.sh

@@ -14,3 +14,7 @@ rm -r data/cutout_cdn/*
 mv data/2010ag/ data/cutout_cdn/
 mv data/2010ai/ data/cutout_cdn/
 mv data/2010h/ data/cutout_cdn/
+
+mv data/sed_output/2010h/ data/
+rm -r data/sed_output/*
+mv data/2010h data/sed_output/

app/host/base_tasks.py

@@ -139,11 +139,11 @@ def find_register_items_meeting_prerequisites(self):
         task_register = TaskRegister.objects.all()
 
         for task_name, status_message in self.prerequisites.items():
-            task = Task.objects.get(name__exact=task_name)
+            task_prereq = Task.objects.get(name__exact=task_name)
             status = Status.objects.get(message__exact=status_message)
 
             current_transients = current_transients & Transient.objects.filter(
-                taskregister__task=task, taskregister__status=status
+                taskregister__task=task_prereq, taskregister__status=status
             )
 
         return task_register.filter(transient__in=list(current_transients), task=task)

app/host/fixtures/example/2010h.yaml

@@ -528,3 +528,43 @@
       flux_error: 2666271.5
       magnitude: 0.0
       magnitude_error: 0.0
+- model: host.SEDFittingResult
+  fields:
+      transient: [2010h]
+      aperture: [2010h_PanSTARRS_g_global]
+      posterior: /sed_output/2010h/2010h_global.h5
+      log_mass_16: 11.11623006468808
+      log_mass_50: 11.119444597406906
+      log_mass_84: 11.122195767038848
+      log_sfr_16: -9.885957158543373
+      log_sfr_50: -9.452641904357886
+      log_sfr_84: -7.854718482318168
+      log_ssfr_16: -21.005705687227554
+      log_ssfr_50: -20.56957875030653
+      log_ssfr_84: -18.975108348918212
+      log_age_16: 3.3558538428424978
+      log_age_50: 3.3848699023305824
+      log_age_84: 3.4566111510622477
+      log_tau_16: 0.10096442863244388
+      log_tau_50: 0.1046197140354783
+      log_tau_84: 0.1261118048050916
+- model: host.SEDFittingResult
+  fields:
+      transient: [2010h]
+      aperture: [2010h_local]
+      posterior: /sed_output/2010h/2010h_local.h5
+      log_mass_16: 8.849994109986605
+      log_mass_50: 8.894097969541852
+      log_mass_84: 8.96004598530449
+      log_sfr_16: -8.953639946834054
+      log_sfr_50: -5.493737812812891
+      log_sfr_84: -2.949760445521605
+      log_ssfr_16: -17.851239889289612
+      log_ssfr_50: -14.381920176924307
+      log_ssfr_84: -11.890886258100974
+      log_age_16: 2.7893040809653034
+      log_age_50: 3.9156213648412592
+      log_age_84: 6.416906088053504
+      log_tau_16: 0.1255658907438466
+      log_tau_50: 0.22027085216263398
+      log_tau_84: 0.5961470357024149

app/host/migrations/0008_auto_20220822_2335.py

@@ -0,0 +1,76 @@
+# Generated by Django 3.2.9 on 2022-08-22 23:35
+import django.db.models.deletion
+import host.models
+from django.db import migrations
+from django.db import models
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ("host", "0007_transient_processing_status"),
+    ]
+
+    operations = [
+        migrations.CreateModel(
+            name="SEDFittingResult",
+            fields=[
+                (
+                    "id",
+                    models.BigAutoField(
+                        auto_created=True,
+                        primary_key=True,
+                        serialize=False,
+                        verbose_name="ID",
+                    ),
+                ),
+                (
+                    "posterior",
+                    models.FileField(
+                        blank=True, null=True, upload_to=host.models.hdf5_file_path
+                    ),
+                ),
+                ("log_mass_16", models.FloatField(blank=True, null=True)),
+                ("log_mass_50", models.FloatField(blank=True, null=True)),
+                ("log_mass_84", models.FloatField(blank=True, null=True)),
+                ("log_sfr_16", models.FloatField(blank=True, null=True)),
+                ("log_sfr_50", models.FloatField(blank=True, null=True)),
+                ("log_sfr_84", models.FloatField(blank=True, null=True)),
+                ("log_ssfr_16", models.FloatField(blank=True, null=True)),
+                ("log_ssfr_50", models.FloatField(blank=True, null=True)),
+                ("log_ssfr_84", models.FloatField(blank=True, null=True)),
+                ("log_age_16", models.FloatField(blank=True, null=True)),
+                ("log_age_50", models.FloatField(blank=True, null=True)),
+                ("log_age_84", models.FloatField(blank=True, null=True)),
+                ("log_tau_16", models.FloatField(blank=True, null=True)),
+                ("log_tau_50", models.FloatField(blank=True, null=True)),
+                ("log_tau_84", models.FloatField(blank=True, null=True)),
+                (
+                    "aperture",
+                    models.ForeignKey(
+                        blank=True,
+                        null=True,
+                        on_delete=django.db.models.deletion.CASCADE,
+                        to="host.aperture",
+                    ),
+                ),
+                (
+                    "transient",
+                    models.ForeignKey(
+                        blank=True,
+                        null=True,
+                        on_delete=django.db.models.deletion.CASCADE,
+                        to="host.transient",
+                    ),
+                ),
+            ],
+        ),
+        migrations.AddField(
+            model_name="host",
+            name="photometric_redshift",
+            field=models.FloatField(blank=True, null=True),
+        ),
+        migrations.DeleteModel(
+            name="ProspectorResult",
+        ),
+    ]

app/host/models.py

@@ -77,6 +77,7 @@ class Host(SkyObject):
 
     name = models.CharField(max_length=100, blank=True, null=True)
     redshift = models.FloatField(null=True, blank=True)
+    photometric_redshift = models.FloatField(null=True, blank=True)
     milkyway_dust_reddening = models.FloatField(null=True, blank=True)
     objects = HostManager()
 
@@ -255,7 +256,9 @@ def transmission_curve(self):
 
         wavelength = transmission_curve[0].to_numpy()
         transmission = transmission_curve[1].to_numpy()
-        return observate.Filter(nick=self.name, data=(wavelength, transmission))
+        return observate.Filter(
+            kname=self.name, nick=self.name, data=(wavelength, transmission)
+        )
 
 
 class Catalog(models.Model):
@@ -281,13 +284,20 @@ class CatalogPhotometry(models.Model):
     catalog = models.ForeignKey(Catalog, on_delete=models.CASCADE)
 
 
-def fits_file_path(instance, filename):
+def fits_file_path(instance):
     """
     Constructs a file path for a fits image
     """
     return f"{instance.host}/{instance.filter.survey}/{instance.filter}.fits"
 
 
+def hdf5_file_path(instance):
+    """
+    Constructs a file path for a HDF5 image
+    """
+    return f"{instance.transient.name}/{instance.transient.name}_{instance.aperture.type}.h5"
+
+
 class Cutout(models.Model):
     """
     Model to represent a cutout image of a host galaxy
@@ -368,19 +378,38 @@ def flux_error_rounded(self):
         return round(self.flux_error, 2)
 
 
-class ProspectorResult(models.Model):
+class SEDFittingResult(models.Model):
     """Model to store prospector results"""
 
-    host = models.ForeignKey(Host, on_delete=models.CASCADE, null=True, blank=True)
-    posterior = models.FileField(upload_to=fits_file_path, null=True, blank=True)
+    transient = models.ForeignKey(
+        Transient, on_delete=models.CASCADE, null=True, blank=True
+    )
+    aperture = models.ForeignKey(
+        Aperture, on_delete=models.CASCADE, null=True, blank=True
+    )
+    posterior = models.FileField(upload_to=hdf5_file_path, null=True, blank=True)
     log_mass_16 = models.FloatField(null=True, blank=True)
     log_mass_50 = models.FloatField(null=True, blank=True)
     log_mass_84 = models.FloatField(null=True, blank=True)
 
+    log_sfr_16 = models.FloatField(null=True, blank=True)
+    log_sfr_50 = models.FloatField(null=True, blank=True)
+    log_sfr_84 = models.FloatField(null=True, blank=True)
+
     log_ssfr_16 = models.FloatField(null=True, blank=True)
     log_ssfr_50 = models.FloatField(null=True, blank=True)
     log_ssfr_84 = models.FloatField(null=True, blank=True)
 
+    # SFR ~ exp(-age/tau)
+    # https://prospect.readthedocs.io/en/latest/sfhs.html?highlight=tau#parametric-sfh
+    log_age_16 = models.FloatField(null=True, blank=True)
+    log_age_50 = models.FloatField(null=True, blank=True)
+    log_age_84 = models.FloatField(null=True, blank=True)
+
+    log_tau_16 = models.FloatField(null=True, blank=True)
+    log_tau_50 = models.FloatField(null=True, blank=True)
+    log_tau_84 = models.FloatField(null=True, blank=True)
+
 
 class TaskRegisterSnapshot(models.Model):
     """

app/host/photometric_calibration.py

@@ -57,5 +57,19 @@ def jansky_to_maggies(flux_density_jansky):
     return flux_density_jansky / 3631.0
 
 
+def mJy_to_maggies(flux_mJy):
+    """
+    Converts spectral flux density from mJy to units of maggies.
+    """
+    return flux_mJy * 10 ** (-0.4 * 23.9)
+
+
+def maggies_to_mJy(flux_maggies):
+    """
+    Converts spectral flux density from maggies to units of mJy.
+    """
+    return flux_maggies * 10 ** (0.4 * 23.9)
+
+
 def magnitude_to_flux_density(magnitude, survey):
     return 0.0

app/host/plotting_utils.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import pandas as pd
+import prospect.io.read_results as reader
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from astropy.visualization import AsinhStretch
@@ -25,6 +26,8 @@
 from bokeh.transform import cumsum
 from host.catalog_photometry import filter_information
 from host.host_utils import survey_list
+from host.photometric_calibration import maggies_to_mJy
+from host.prospector import build_obs  # , build_model
 
 from .models import Aperture
 
@@ -184,7 +187,7 @@ def plot_cutout_image(
     return {"bokeh_cutout_script": script, "bokeh_cutout_div": div}
 
 
-def plot_sed(aperture_photometry=None, type=""):
+def plot_sed(aperture_photometry=None, sed_results_file=None, type=""):
     """
     Plot SED from aperture photometry.
     """
@@ -207,14 +210,25 @@ def plot_sed(aperture_photometry=None, type=""):
         width=700,
         height=400,
         min_border=0,
-        toolbar_location=None,
+        #    toolbar_location=None,
         x_axis_type="log",
         x_axis_label="Wavelength [Angstrom]",
         y_axis_label="Flux",
     )
 
     fig = plot_errorbar(fig, wavelength, flux, yerr=flux_error)
 
+    if sed_results_file is not None:
+        print(sed_results_file)
+        result, obs, model = reader.results_from(sed_results_file, dangerous=False)
+
+        best = result["bestfit"]
+        a = result["obs"]["redshift"] + 1
+        fig.line(a * best["restframe_wavelengths"], maggies_to_mJy(best["spectrum"]))
+        if obs["filters"] is not None:
+            pwave = [f.wave_effective for f in obs["filters"]]
+            fig.circle(pwave, maggies_to_mJy(best["photometry"]))
+
     # xaxis = LinearAxis()
     # figure.add_layout(xaxis, 'below')