Merge pull request #16 from NeoGeographyToolkit/estimate_age_and_crat…

…er_d_fix

Crater diameter and estimate_age fix

src/python/synthterrain/crater/__init__.py

@@ -42,12 +42,7 @@ def synthesize(
        synthesized from the input parameters.
     """
     if production_fn is None:
-        if crater_dist.a >= 10:
-            production_fn = functions.NPF(a=crater_dist.a, b=crater_dist.b)
-        elif crater_dist.b <= 2.76:
-            production_fn = functions.Grun(a=crater_dist.a, b=crater_dist.b)
-        else:
-            production_fn = functions.GNPF(a=crater_dist.a, b=crater_dist.b)
+        production_fn = determine_production_function(crater_dist.a, crater_dist.b)
     logger.info(f"Production function is {production_fn.__class__}")
 
     # Get craters
@@ -90,6 +85,15 @@ def synthesize(
     return df
 
 
+def determine_production_function(a: float, b: float):
+    if a >= 10:
+        return functions.NPF(a, b)
+    elif b <= 2.76:
+        return functions.Grun(a, b)
+    else:
+        return functions.GNPF(a, b)
+
+
 def random_points(poly: Polygon, num_points: int):
     """Returns two lists, the first being the x coordinates, and the second
        being the y coordinates, each *num_points* long that represent

src/python/synthterrain/crater/age.py

@@ -70,9 +70,9 @@ def estimate_age(diameter, dd, max_age):
 
 def estimate_age_by_bin(
     df,
-    pd_csfd,
-    eq_csfd,
     num=50,
+    pd_csfd=None,
+    eq_csfd=None,
 ) -> pd.DataFrame:
     """
     Returns a pandas DataFrame identical to the input *df* but with the
@@ -85,50 +85,66 @@ def estimate_age_by_bin(
     scale boundaries (using the numpy geomspace() function) between the maximum
     and minimum diameter values.
 
-    Then, the center diameter of each bin has its equilibrium age determined via
-    equilibrium_ages() and diffuse_d_over_D() is run to evaluate the d/D ratio
-    of that crater over its lifetime.  Then, for each crater in the bin, its d/D
-    is compared to the d/D ratios from the diffusion run, and an estimated age
-    is assigned.
+    Then, the center diameter of each bin has diffuse_d_over_D()
+    run to evaluate the d/D ratio of that crater over the lifetime
+    of the solar system (if *pd_csfd* and *eq_csfd* are both None), or
+    the diffusion calculation is run for the appropriate equilibrium age
+    of each crater diameter.
+
+    Then, for each crater in the bin, its d/D is compared to the d/D ratios from the
+    diffusion run, and an estimated age is assigned.
     """
     if df.empty:
         raise ValueError("The provided dataframe has no rows.")
 
     logger.info("estimate_age_by_bin start.")
-    bin_edges = np.geomspace(
-        df["diameter"].min(), df["diameter"].max(), num=num + 1
-    )
-    # logger.info(f"{df.shape[0]} craters")
-    logger.info(
-        f"Divided the craters into {num} diameter bins (not all bins may have "
-        "craters)"
+    if df["diameter"].min() == df["diameter"].max():
+        bin_edges = 1
+        total_bins = 1
+        logger.info("The craters are all the same size.")
+    else:
+        bin_edges = np.geomspace(
+            df["diameter"].min(), df["diameter"].max(), num=num + 1
+        )
+        total_bins = num
+        # logger.info(f"{df.shape[0]} craters")
+        logger.info(
+            f"Divided the craters into {num} diameter bins (not all bins may have "
+            "craters)"
+        )
+
+    df["diameter_bin"] = pd.cut(
+        df["diameter"], bins=bin_edges, include_lowest=True,
     )
-    df["diameter_bin"] = pd.cut(df["diameter"], bins=bin_edges, include_lowest=True)
+
     # df["equilibrium_age"] = equilibrium_ages(df["diameter"], pd_csfd, eq_csfd)
     df["age"] = 0
 
     for i, (interval, count) in enumerate(
             df["diameter_bin"].value_counts(sort=False).items()
     ):
         logger.info(
-            f"Processing bin {i}/{num}, interval: {interval}, count: {count}"
+            f"Processing bin {i + 1}/{total_bins}, interval: {interval}, count: {count}"
         )
 
         if count == 0:
             continue
 
         fresh_dd = stopar_fresh_dd(interval.mid)
 
-        eq_age = equilibrium_age(interval.mid, pd_csfd, eq_csfd)
+        if pd_csfd is not None and eq_csfd is not None:
+            age = equilibrium_age(interval.mid, pd_csfd, eq_csfd)
+        else:
+            age = 4.5e9
 
         dd_rev_list = list(reversed(diffuse_d_over_D(
             interval.mid,
-            eq_age,
+            age,
             start_dd_adjust=fresh_dd,
             return_steps=True
         )))
         nsteps = len(dd_rev_list)
-        years_per_step = eq_age / nsteps
+        years_per_step = age / nsteps
 
         def guess_age(dd):
             age_step = bisect.bisect_left(dd_rev_list, dd)

src/python/synthterrain/crater/cli_convert.py

@@ -40,6 +40,15 @@ def arg_parser():
              "craters may still yield a zero age if the d/D ratio was large relative "
              "to the diameter, indicating a very fresh crater."
     )
+    parser.add_argument(
+        "--full_age",
+        action="store_true",
+        help="Ignored unless --estimate_ages is also given.  When provided, it will "
+             "cause the diffusion calculation to run for the age of the solar system "
+             "instead of just the equilibrium age for each crater size.  This may "
+             "provide improved age estimates, but could also cause longer run times. "
+             "Please use with caution."
+    )
     parser.add_argument(
         "infile",
         type=Path,
@@ -77,17 +86,22 @@ def main():
         )
 
     if args.estimate_ages:
-        a = df["diameter"].min()
-        b = df["diameter"].max()
-        pd_func = crater_func.GNPF(a=a, b=b)
-        eq_func = crater_func.VIPER_Env_Spec(a=a, b=b)
+        if args.full_age:
+            pd_csfd = None
+            eq_csfd = None
+        else:
+            a = df["diameter"].min()
+            b = df["diameter"].max()
+            pd_csfd = crater.determine_production_function(a=a, b=b).csfd
+            eq_csfd = crater_func.VIPER_Env_Spec(a=a, b=b).csfd
+
         try:
             if "age" in df.columns:
                 df[df["age"] == 0] = estimate_age_by_bin(
-                    df[df["age"] == 0], pd_func.csfd, eq_func.csfd, num=50
+                    df[df["age"] == 0], 50, pd_csfd, eq_csfd
                 )
             else:
-                df = estimate_age_by_bin(df, pd_func.csfd, eq_func.csfd, num=50)
+                df = estimate_age_by_bin(df, 50, pd_csfd, eq_csfd)
         except ValueError:
             logger.error(
                 "The provided file has no craters with an age of zero."

tests/python/crater/test_age.py

@@ -50,9 +50,9 @@ def test_estimate_age_by_bin(self):
         )
         df_out = age.estimate_age_by_bin(
             df,
+            50,  # the bin size can have a real impact here.
             pd_func.csfd,
             eq_func.csfd,
-            num=50,  # the bin size can have a real impact here.
         )
 
         age_series = pd.Series([
@@ -61,3 +61,18 @@ def test_estimate_age_by_bin(self):
             8072614], name="age")
 
         pd.testing.assert_series_equal(age_series, df_out["age"])
+
+
+        df2 = pd.DataFrame(data={
+            'diameter': [100., 100., 100., 100.],
+            'd/D': [0.01, 0.06, 0.10, 0.17]
+        })
+
+        df2_out = age.estimate_age_by_bin(
+            df2,
+            50,  # With only one diameter, num is irrelevant
+        )
+
+        age2_series = pd.Series([4500000000, 4500000000, 1388888888, 0], name="age")
+
+        pd.testing.assert_series_equal(age2_series, df2_out["age"])