diff --git a/snsim/geo_utils.py b/snsim/geo_utils.py index 10c3215..1249695 100644 --- a/snsim/geo_utils.py +++ b/snsim/geo_utils.py @@ -17,6 +17,7 @@ def _format_corner(corner, RA): # - RA_0 < RA_1 # - RA_3 < RA_2 # - RA_0 and RA_3 on the same side of the field center + # corner[fields, corner, subfields, ra/dec] sign = (corner[:, 3, :, 0] - RA[:, None]) * (corner[:, 0, :, 0] - RA[:, None]) < 0 comp = corner[:, 0, :, 0] < corner[:, 3, :, 0] @@ -46,7 +47,12 @@ def _compute_area(polygon): def _compute_polygon(corners): - """Create polygon on a sphere, check for edges conditions.""" + """Create polygon on a sphere, check for edges conditions. + + Notes + ----- + corners[corner, subfields, ra/dec] + """ # Create polygons polygons = gpd.GeoSeries([shp_geo.Polygon(corners[:, j, :]) for j in range(corners.shape[1])]) diff --git a/snsim/tests/test_geout.py b/snsim/tests/test_geout.py new file mode 100644 index 0000000..ae88de6 --- /dev/null +++ b/snsim/tests/test_geout.py @@ -0,0 +1,26 @@ +import numpy as np +import shapely.geometry as shp_geo +from numpy.testing import assert_almost_equal +from snsim import geo_utils as geo_ut + +def test_compute_area(): + corners = np.array([[[0, np.pi / 2]], [[2 * np.pi, np.pi / 2]] , + [[2 * np.pi, -np.pi / 2]], [[0, -np.pi / 2]]]) + + polygon = geo_ut._compute_polygon(corners) + + area = geo_ut._compute_area(polygon) + assert_almost_equal(area, 4 * np.pi) + + +def test_compute_polygon(): + + # Regular polygon + corners = np.array([[[np.pi - 0.02, 0.02]], [[np.pi + 0.02, 0.02]] , + [[np.pi + 0.02, -0.02]], [[np.pi - 0.02, -0.02]]]) + + polygon = geo_ut._compute_polygon(corners) + + ra = np.random.uniform(0, 2 * np.pi, 10 + + \ No newline at end of file