Skip to content

Commit ecc6a27

Browse files
olivialynnolivialynn
committed
Split read/write/verify into 2 different notebooks
1 parent b039bbe commit ecc6a27

2 files changed

Lines changed: 352 additions & 279 deletions

File tree

src/skymap_to_poly_coords/01 - Storing skymaps.ipynb renamed to src/skymap_to_poly_coords/01 - Convert and write poly skymap.ipynb

Lines changed: 5 additions & 279 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
},
1919
{
2020
"cell_type": "code",
21-
"execution_count": 1,
21+
"execution_count": null,
2222
"id": "b046e4b1",
2323
"metadata": {},
2424
"outputs": [],
@@ -27,13 +27,8 @@
2727
"from pathlib import Path\n",
2828
"import math\n",
2929
"import numpy as np\n",
30-
"\n",
31-
"\n",
32-
"from lsst.sphgeom import Box, ConvexPolygon, UnitVector3d\n",
33-
"from math import pi\n",
34-
"from lsst.sphgeom import LonLat, UnitVector3d, ConvexPolygon\n",
35-
"import yaml\n",
36-
"from lsst.sphgeom import UnitVector3d, ConvexPolygon"
30+
"from lsst.sphgeom import Box, ConvexPolygon, LonLat, UnitVector3d\n",
31+
"import yaml"
3732
]
3833
},
3934
{
@@ -63,171 +58,13 @@
6358
"lsst_skymap"
6459
]
6560
},
66-
{
67-
"cell_type": "markdown",
68-
"id": "c035d29d",
69-
"metadata": {},
70-
"source": [
71-
"## Utility functions for later"
72-
]
73-
},
74-
{
75-
"cell_type": "markdown",
76-
"id": "95b64469",
77-
"metadata": {},
78-
"source": [
79-
"### Check if a point is in a tract (using lsst.skymap)"
80-
]
81-
},
82-
{
83-
"cell_type": "code",
84-
"execution_count": 3,
85-
"id": "2f076531",
86-
"metadata": {},
87-
"outputs": [],
88-
"source": [
89-
"def box_to_convex_polygon(box: Box) -> ConvexPolygon:\n",
90-
" if box.isEmpty():\n",
91-
" raise ValueError(\"Cannot convert an empty Box to a ConvexPolygon.\")\n",
92-
"\n",
93-
" # Get the corners of the box\n",
94-
" lon_a, lon_b = box.getLon().getA().asRadians(), box.getLon().getB().asRadians()\n",
95-
" lon_min = min(lon_a, lon_b)\n",
96-
" lon_max = max(lon_a, lon_b)\n",
97-
" lat_a, lat_b = box.getLat().getA().asRadians(), box.getLat().getB().asRadians()\n",
98-
" lat_min = min(lat_a, lat_b)\n",
99-
" lat_max = max(lat_a, lat_b)\n",
100-
" # todo : this may be an improper assumption, considering RA wrap around!!\n",
101-
"\n",
102-
" bottom_left = LonLat.fromRadians(lon_min, lat_min)\n",
103-
" bottom_right = LonLat.fromRadians(lon_max, lat_min)\n",
104-
" top_right = LonLat.fromRadians(lon_max, lat_max)\n",
105-
" top_left = LonLat.fromRadians(lon_min, lat_max)\n",
106-
"\n",
107-
" # Convert corners to UnitVector3d\n",
108-
" vertices = [\n",
109-
" UnitVector3d(bottom_left),\n",
110-
" UnitVector3d(bottom_right),\n",
111-
" UnitVector3d(top_right),\n",
112-
" UnitVector3d(top_left),\n",
113-
" ]\n",
114-
"\n",
115-
" # Create and return the ConvexPolygon\n",
116-
" return ConvexPolygon(vertices)"
117-
]
118-
},
119-
{
120-
"cell_type": "markdown",
121-
"id": "7ce7f70c",
122-
"metadata": {},
123-
"source": [
124-
"### Get a poly from a tract ID"
125-
]
126-
},
127-
{
128-
"cell_type": "code",
129-
"execution_count": 4,
130-
"id": "6eb59d04",
131-
"metadata": {},
132-
"outputs": [],
133-
"source": [
134-
"def get_poly_from_tract_id(tract_id, inner=False) -> ConvexPolygon:\n",
135-
" tract = lsst_skymap.generateTract(tract_id)\n",
136-
" if inner:\n",
137-
" res = tract.inner_sky_region\n",
138-
" else:\n",
139-
" res = tract.outer_sky_polygon\n",
140-
" if isinstance(res, Box):\n",
141-
" res = box_to_convex_polygon(res)\n",
142-
" return res"
143-
]
144-
},
145-
{
146-
"cell_type": "markdown",
147-
"id": "c85b5340",
148-
"metadata": {},
149-
"source": [
150-
"### Check if a point (in ra/dec, degrees) is in a given poly"
151-
]
152-
},
153-
{
154-
"cell_type": "code",
155-
"execution_count": 5,
156-
"id": "ce01c001",
157-
"metadata": {},
158-
"outputs": [
159-
{
160-
"data": {
161-
"text/plain": [
162-
"True"
163-
]
164-
},
165-
"execution_count": 5,
166-
"metadata": {},
167-
"output_type": "execute_result"
168-
}
169-
],
170-
"source": [
171-
"\n",
172-
"def point_in_poly(polygon, ra_degrees, dec_degrees):\n",
173-
" vec = UnitVector3d(LonLat.fromDegrees(ra_degrees, dec_degrees))\n",
174-
" return polygon.contains(vec)\n",
175-
"\n",
176-
"\n",
177-
" \n",
178-
"#check_point(1, 0.0, -88.0, inner=False)\n",
179-
"poly = get_poly_from_tract_id(1)\n",
180-
"point_in_poly(poly, 0.0, -88.0)"
181-
]
182-
},
183-
{
184-
"cell_type": "markdown",
185-
"id": "2b30451b",
186-
"metadata": {},
187-
"source": [
188-
"### Check if two polys are equivalent (within tolerance)"
189-
]
190-
},
191-
{
192-
"cell_type": "code",
193-
"execution_count": 6,
194-
"id": "691ac997",
195-
"metadata": {},
196-
"outputs": [],
197-
"source": [
198-
"\n",
199-
"def polys_are_equiv(poly_a, poly_b, rtol=1e-12, atol=1e-14):\n",
200-
" \"\"\"Check if two ConvexPolygons are equivalent within floating point tolerance.\n",
201-
"\n",
202-
" Parameters\n",
203-
" ----------\n",
204-
" poly_a, poly_b : sphgeom.ConvexPolygon\n",
205-
" The polygons to compare.\n",
206-
" rtol : float\n",
207-
" Relative tolerance for np.allclose.\n",
208-
" atol : float\n",
209-
" Absolute tolerance for np.allclose.\n",
210-
"\n",
211-
" Returns\n",
212-
" -------\n",
213-
" bool\n",
214-
" True if all vertices match within tolerance.\n",
215-
" \"\"\"\n",
216-
" verts_a = poly_a.getVertices()\n",
217-
" verts_b = poly_b.getVertices()\n",
218-
"\n",
219-
" if len(verts_a) != len(verts_b):\n",
220-
" return False\n",
221-
"\n",
222-
" return np.allclose(verts_a, verts_b, rtol=rtol, atol=atol)"
223-
]
224-
},
22561
{
22662
"cell_type": "markdown",
22763
"id": "3ce11262",
22864
"metadata": {},
22965
"source": [
230-
"*A note on rings sky map pixelization (in html comment in this cell)*\n",
66+
"## A note on rings sky map pixelization \n",
67+
"*(in html comment in this cell)*\n",
23168
"<!---## Rings sky map pixelization\n",
23269
"*From the RingsSkyMap docstring in lsst.skymap:*\n",
23370
"\n",
@@ -387,117 +224,6 @@
387224
"write_polygons(lsst_skymap, outer_poly_path, inner=False)"
388225
]
389226
},
390-
{
391-
"cell_type": "markdown",
392-
"id": "4e89bd8f",
393-
"metadata": {},
394-
"source": [
395-
"### Read inner_poly and outer_poly"
396-
]
397-
},
398-
{
399-
"cell_type": "code",
400-
"execution_count": 11,
401-
"id": "7c0022ee",
402-
"metadata": {},
403-
"outputs": [],
404-
"source": [
405-
"def load_polygons(yaml_path):\n",
406-
" \"\"\"Load exact inner or outer polygons from a YAML file using 3D unit vectors.\n",
407-
"\n",
408-
" Parameters\n",
409-
" ----------\n",
410-
" yaml_path : str\n",
411-
" Path to the YAML file written by `write_polygons`.\n",
412-
"\n",
413-
" Returns\n",
414-
" -------\n",
415-
" dict\n",
416-
" Mapping from tract ID (int) to sphgeom.ConvexPolygon.\n",
417-
" \"\"\"\n",
418-
" with open(yaml_path, \"r\") as f:\n",
419-
" data = yaml.safe_load(f)\n",
420-
"\n",
421-
" poly_dict = {}\n",
422-
"\n",
423-
" for tract_id_str, vec_list in data[\"tracts\"].items():\n",
424-
" tract_id = int(tract_id_str)\n",
425-
"\n",
426-
" unit_vecs = [UnitVector3d(*vec) for vec in vec_list]\n",
427-
"\n",
428-
" # Skip degenerate polygons (fewer than 3 unique vertices)\n",
429-
" unique_vecs = {tuple(round(x, 12) for x in v) for v in unit_vecs}\n",
430-
" if len(unique_vecs) < 3:\n",
431-
" print(f\"⚠️ Skipping degenerate tract {tract_id}\")\n",
432-
" continue\n",
433-
"\n",
434-
" poly = ConvexPolygon(unit_vecs)\n",
435-
" poly_dict[tract_id] = poly\n",
436-
"\n",
437-
" print(f\"✅ Loaded {len(poly_dict)} polygons from {yaml_path}\")\n",
438-
" return poly_dict"
439-
]
440-
},
441-
{
442-
"cell_type": "code",
443-
"execution_count": 12,
444-
"id": "3f7e50e9",
445-
"metadata": {},
446-
"outputs": [
447-
{
448-
"name": "stdout",
449-
"output_type": "stream",
450-
"text": [
451-
"⚠️ Skipping degenerate tract 0\n",
452-
"⚠️ Skipping degenerate tract 18937\n",
453-
"✅ Loaded 18936 polygons from /sdf/home/o/olynn/skymap-to-poly-coords/skymaps_out/inner_polys.yaml\n",
454-
"✅ Loaded 18938 polygons from /sdf/home/o/olynn/skymap-to-poly-coords/skymaps_out/outer_polys.yaml\n"
455-
]
456-
}
457-
],
458-
"source": [
459-
"inner_poly_map = load_polygons(inner_poly_path)\n",
460-
"outer_poly_map = load_polygons(outer_poly_path)"
461-
]
462-
},
463-
{
464-
"cell_type": "code",
465-
"execution_count": 16,
466-
"id": "7534eedf",
467-
"metadata": {},
468-
"outputs": [],
469-
"source": [
470-
"# Just patch the polar caps in for now (todo)\n",
471-
"\n",
472-
"inner_poly_map[0] = get_poly_from_tract_id(0, inner=True)\n",
473-
"inner_poly_map[len(inner_poly_map)-1] = get_poly_from_tract_id(len(inner_poly_map)-1, inner=True)"
474-
]
475-
},
476-
{
477-
"cell_type": "markdown",
478-
"id": "4d1fe8bb",
479-
"metadata": {},
480-
"source": [
481-
"### Check our saved-and-loaded tracts against the tracts we read via the lsst.skymap package"
482-
]
483-
},
484-
{
485-
"cell_type": "code",
486-
"execution_count": 17,
487-
"id": "79f1b755",
488-
"metadata": {},
489-
"outputs": [],
490-
"source": [
491-
"tracts_to_check = np.linspace(0, lsst_skymap._numTracts-1, 1000, dtype=int)\n",
492-
"for tract_id in tracts_to_check:\n",
493-
" ground_truth_poly = get_poly_from_tract_id(tract_id, inner=True)\n",
494-
" loaded_poly = inner_poly_map[tract_id]\n",
495-
" if not polys_are_equiv(ground_truth_poly, loaded_poly):\n",
496-
" print(f\"Tract {tract_id} polygons are NOT equivalent!\")\n",
497-
" else:\n",
498-
" continue"
499-
]
500-
},
501227
{
502228
"cell_type": "markdown",
503229
"id": "eca14065",

0 commit comments

Comments
 (0)