Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update convex hull code #789

Closed
wants to merge 12 commits into from
Prev Previous commit
Next Next commit
codestyle check 4
  • Loading branch information
christinawlindberg committed May 26, 2023
commit fe8373a48afc127ae1a923f418c25ed8554fb595
5 changes: 4 additions & 1 deletion beast/tools/cut_catalogs.py
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@
# do convex hull
ra_col = [x for x in phot_cat.colnames if "RA" in x.upper()][0]
dec_col = [x for x in phot_cat.colnames if "DEC" in x.upper()][0]
phot_cat_boundary = convexhull_path(new_phot_cat[ra_col], new_phot_cat[dec_col],

Check warning on line 110 in beast/tools/cut_catalogs.py

View check run for this annotation

Codecov / codecov/patch

beast/tools/cut_catalogs.py#L110

Added line #L110 was not covered by tests
concave_alpha, input_ast_file)
# check if points are inside
ra_col = [x for x in ast_cat.colnames if "RA" in x.upper()][0]
Expand Down Expand Up @@ -258,18 +258,21 @@
[x_coord, y_coord]
).T # there's a weird astropy datatype issue that requires numpy coercion

print("using concave hull, alpha: " + str(concave_alpha) + " (default: 100)")
concave_hull = alpha_shape(coords, concave_alpha)
path_object = Path(np.array(concave_hull.boundary.xy).T)

Check warning on line 263 in beast/tools/cut_catalogs.py

View check run for this annotation

Codecov / codecov/patch

beast/tools/cut_catalogs.py#L261-L263

Added lines #L261 - L263 were not covered by tests

# create diagnostic plot
ax = plt.subplots(1, 1, figsize=(5, 4))
fig = plt.figure(figsize=(5, 4))
ax = fig.subplots(1, 1)
ax.scatter(x_coord, y_coord, s=1)
ax.plot(path_object.vertices.T[0], path_object.vertices.T[1], c="r")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title(r"$\alpha$ = " + str(concave_alpha))
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.savefig(str(input_ast_file) + "_plot.png", bbox_inches="tight")

Check warning on line 275 in beast/tools/cut_catalogs.py

View check run for this annotation

Codecov / codecov/patch

beast/tools/cut_catalogs.py#L266-L275

Added lines #L266 - L275 were not covered by tests

return path_object

Expand All @@ -284,28 +287,28 @@
don't fall inward as much as larger numbers.
Too large, and you lose everything!
"""
if len(coords) < 4:

Check warning on line 290 in beast/tools/cut_catalogs.py

View check run for this annotation

Codecov / codecov/patch

beast/tools/cut_catalogs.py#L290

Added line #L290 was not covered by tests
# When you have a triangle, there is no sense
# in computing an alpha shape.
return geometry.MultiPoint(list(coords)).convex_hull

Check warning on line 293 in beast/tools/cut_catalogs.py

View check run for this annotation

Codecov / codecov/patch

beast/tools/cut_catalogs.py#L293

Added line #L293 was not covered by tests

tri = Delaunay(coords)
triangles = coords[tri.simplices]
a = ((triangles[:, 0, 0] - triangles[:, 1, 0]) ** 2 + (triangles[:, 0, 1] - triangles[:, 1, 1]) ** 2) ** 0.5
b = ((triangles[:, 1, 0] - triangles[:, 2, 0]) ** 2 + (triangles[:, 1, 1] - triangles[:, 2, 1]) ** 2) ** 0.5
c = ((triangles[:, 2, 0] - triangles[:, 0, 0]) ** 2 + (triangles[:, 2, 1] - triangles[:, 0, 1]) ** 2) ** 0.5
s = (a + b + c) / 2.0
areas = (s * (s - a) * (s - b) * (s - c)) ** 0.5
circums = a * b * c / (4.0 * areas)
filtered = triangles[circums < (1.0 / alpha)]
edge1 = filtered[:, (0, 1)]
edge2 = filtered[:, (1, 2)]
edge3 = filtered[:, (2, 0)]
edge_points = np.unique(np.concatenate((edge1, edge2, edge3)), axis=0).tolist()
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))

Check warning on line 309 in beast/tools/cut_catalogs.py

View check run for this annotation

Codecov / codecov/patch

beast/tools/cut_catalogs.py#L295-L309

Added lines #L295 - L309 were not covered by tests

return unary_union(triangles)

Check warning on line 311 in beast/tools/cut_catalogs.py

View check run for this annotation

Codecov / codecov/patch

beast/tools/cut_catalogs.py#L311

Added line #L311 was not covered by tests


if __name__ == "__main__": # pragma: no cover
Expand Down