Skip to content

Commit

Permalink
Merge pull request #244 from astronomy-commons/sean/margin-notebook
Browse files Browse the repository at this point in the history 
Add margin documentation notebook
  • Loading branch information
@smcguire-cmu
smcguire-cmu authored Apr 2, 2024
2 parents 1584cff + bda5722 commit ebcabd1
Show file tree
Hide file tree
Showing 3 changed files with 377 additions and 3 deletions.
Binary file added docs/tutorials/_static/margin-pix.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/tutorials/_static/pixel-boundary-example.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
380 changes: 377 additions & 3 deletions docs/tutorials/margins.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,16 +4,390 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"# Margins"
"# Margins\n",
"\n",
"LSDB can handle datasets larger than memory by breaking them down into smaller spatially-connected parts and working on each part one at a time. One of the main tasks enabled by LSDB are spatial queries such as cross-matching; To ensure accurate comparisons, all nearby data points need to be loaded simultaneously. LSDB uses HiPSCat's method of organizing data spatially to achieve this. However, there's a limitation: at the boundaries of each divided section, some data points are going to be missed. This means that for operations requiring comparisons with neighboring points, such as cross-matching, the process might miss some matches for points near these boundaries because not all nearby points are included when analyzing one section at a time.\n",
"\n",
"![Margin Boundary Example](_static/pixel-boundary-example.png)\n",
"*Here we see an example of a boundary between HEALPix pixels, where the green points are in one partition and the red points in another. Working with one partition at a time, we would miss potential matches with points close to the boundary*\n",
"\n",
"To solve this, we could try to also load the neighboring partitions for each partition we crossmatch. However, this would mean needing to load lots of unnecessary data, slowing down operations and causing issues with running out of memory. So for each catalog we also create a margin cache. This means that for each partition, we create a file that contains the points in the catalog within a certain distance to the pixel's boundary.\n",
"\n",
"![Margin Cache Example](_static/margin-pix.png)\n",
"*An example of a margin cache (orange) for the same green pixel. The margin cache for this pixel contains the points within 10 arcseconds of the boundary.*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
"metadata": {
"collapsed": false
},
"source": [
"## Loading a Margin Catalog\n",
"\n",
"The margin cache is stored as a separate HiPSCat catalog with the same partitioning as the main catalog.\n",
"\n",
"To load a margin cache, we first load the margin catalog the same as other HiPSCat catalogs with a `lsdb.read_hipscat` call."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-26T20:32:10.418696Z",
"start_time": "2024-03-26T20:32:08.688730Z"
}
},
"outputs": [],
"source": [
"!pip install aiohttp --quiet # for reading from http"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:38:39.721106Z",
"start_time": "2024-03-28T20:38:35.711117Z"
}
},
"outputs": [],
"source": [
"import lsdb\n",
"\n",
"surveys_path = \"https://epyc.astro.washington.edu/~lincc-frameworks/half_degree_surveys\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:39:05.137393Z",
"start_time": "2024-03-28T20:39:02.808295Z"
}
},
"outputs": [],
"source": [
"ztf_margin_path = f\"{surveys_path}/ztf/ztf_object_margin\"\n",
"ztf_margin = lsdb.read_hipscat(ztf_margin_path)\n",
"print(f\"Margin size: {ztf_margin.hc_structure.catalog_info.margin_threshold} arcsec\")\n",
"ztf_margin"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Here we see the margin catalog that has been generated, in this case using a margin threshold of 10 arcseconds."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Then we load the object catalog, setting the `margin_cache` parameter to the margin catalog we have just loaded."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:40:57.666217Z",
"start_time": "2024-03-28T20:40:55.333646Z"
}
},
"outputs": [],
"source": [
"ztf_object_path = f\"{surveys_path}/ztf/ztf_object\"\n",
"ztf_object = lsdb.read_hipscat(ztf_object_path, margin_cache=ztf_margin)\n",
"ztf_object"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"The `ztf_object` catalog has now been loaded with its margin cache (This dataset is a small sample of the full DR14 Catalog). We can plot the catalog and its margin to see this."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:49:47.638073Z",
"start_time": "2024-03-28T20:49:47.616533Z"
}
},
"outputs": [],
"source": [
"# Defining a function to plot the points in a pixel and the pixel boundary\n",
"\n",
"import numpy as np\n",
"from matplotlib.patches import Polygon\n",
"from matplotlib import pyplot as plt\n",
"import healpy as hp\n",
"\n",
"\n",
"def plot_points(\n",
" pixel_dfs, order, pixel, colors, ra_columns, dec_columns, xlim=None, ylim=None, markers=None, alpha=1\n",
"):\n",
" ax = plt.subplot()\n",
"\n",
" # Plot hp pixel bounds\n",
" nsides = hp.order2nside(order)\n",
" pix0_bounds = hp.vec2dir(hp.boundaries(nsides, pixel, step=100, nest=True), lonlat=True)\n",
" lon = pix0_bounds[0]\n",
" lat = pix0_bounds[1]\n",
" vertices = np.vstack([lon.ravel(), lat.ravel()]).transpose()\n",
" p = Polygon(vertices, closed=True, edgecolor=\"#3b81db\", facecolor=\"none\")\n",
" ax.add_patch(p)\n",
"\n",
" if markers is None:\n",
" markers = [\"+\"] * len(pixel_dfs)\n",
"\n",
" # plot the points\n",
" for pixel_df, color, ra_column, dec_column, marker in zip(\n",
" pixel_dfs, colors, ra_columns, dec_columns, markers\n",
" ):\n",
" ax.scatter(\n",
" pixel_df[ra_column].values,\n",
" pixel_df[dec_column].values,\n",
" c=color,\n",
" marker=marker,\n",
" linewidths=1,\n",
" alpha=alpha,\n",
" )\n",
"\n",
" # plotting configuration\n",
" VIEW_MARGIN = 2\n",
" xlim_low = np.min(lon) - VIEW_MARGIN if xlim is None else xlim[0]\n",
" xlim_high = np.max(lon) + VIEW_MARGIN if xlim is None else xlim[1]\n",
" ylim_low = np.min(lat) - VIEW_MARGIN if ylim is None else ylim[0]\n",
" ylim_high = np.max(lat) + VIEW_MARGIN if ylim is None else ylim[1]\n",
"\n",
" plt.xlim(xlim_low, xlim_high)\n",
" plt.ylim(ylim_low, ylim_high)\n",
" plt.xlabel(\"ra\")\n",
" plt.ylabel(\"dec\")\n",
" plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:49:50.702578Z",
"start_time": "2024-03-28T20:49:49.245902Z"
}
},
"outputs": [],
"source": [
"# the healpix pixel to plot\n",
"order = 3\n",
"pixel = 434\n",
"\n",
"# Plot the points from the specified ztf pixel in green, and from the pixel's margin cache in red\n",
"plot_points(\n",
" [\n",
" ztf_object.get_partition(order, pixel).compute(),\n",
" ztf_object.margin.get_partition(order, pixel).compute(),\n",
" ],\n",
" order,\n",
" pixel,\n",
" [\"green\", \"red\"],\n",
" [\"ra\", \"ra\"],\n",
" [\"dec\", \"dec\"],\n",
" xlim=[179.5, 180.1],\n",
" ylim=[9.4, 10.0],\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## Using the Margin Catalog\n",
"\n",
"Performing operations like cross-matching and joining require a margin to be loaded in the catalog on the right side of the operation. If this right catalog has been loaded with a margin, the function will be carried out accurately using the margin, and by default will throw an error if the margin has not been set. This can be overwritten using the `require_right_margin` parameter, but this may cause inaccurate results!\n",
"\n",
"We can see this trying to perform a crossmatch with gaia"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:49:55.217215Z",
"start_time": "2024-03-28T20:49:52.586757Z"
}
},
"outputs": [],
"source": [
"gaia = lsdb.read_hipscat(f\"{surveys_path}/gaia\")\n",
"gaia"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"If we perform a crossmatch with gaia on the left and the ztf catalog we loaded with a margin on the right, the function works and we get the result"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:49:55.477297Z",
"start_time": "2024-03-28T20:49:55.222520Z"
}
},
"outputs": [],
"source": [
"gaia.crossmatch(ztf_object)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"If we try the other way around, we have not loaded the right catalog (gaia) with a margin cache, and so we get an error"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:49:57.849868Z",
"start_time": "2024-03-28T20:49:57.843828Z"
}
},
"outputs": [],
"source": [
"try:\n",
" ztf_object.crossmatch(gaia)\n",
"except ValueError as e:\n",
" print(f\"Error: {e}\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"We can plot the result of the crossmatch below, with the gaia objects in green and the ztf objects in red"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:50:01.240040Z",
"start_time": "2024-03-28T20:49:59.007988Z"
}
},
"outputs": [],
"source": [
"crossmatch_result = gaia.crossmatch(ztf_object)\n",
"\n",
"plot_points(\n",
" [\n",
" crossmatch_result.get_partition(order, pixel).compute(),\n",
" crossmatch_result.get_partition(order, pixel).compute(),\n",
" ],\n",
" order,\n",
" pixel,\n",
" [\"green\", \"red\"],\n",
" [\"ra_gaia_halfdegree\", \"ra_ztf_object_halfdegree\"],\n",
" [\"dec_gaia_halfdegree\", \"dec_ztf_object_halfdegree\"],\n",
" xlim=[179.5, 180.1],\n",
" ylim=[9.4, 10.0],\n",
" markers=[\"+\", \"x\"],\n",
" alpha=0.8,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## Filtering Catalogs\n",
"\n",
"Any filtering operations applied to the catalog are also performed to the margin."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2024-03-28T20:50:05.236767Z",
"start_time": "2024-03-28T20:50:04.070543Z"
}
},
"outputs": [],
"source": [
"small_sky_box_filter = ztf_object.box(ra=[179.9, 180], dec=[9.5, 9.7])\n",
"\n",
"# Plot the points from the filtered ztf pixel in green, and from the pixel's margin cache in red\n",
"plot_points(\n",
" [\n",
" small_sky_box_filter.get_partition(order, pixel).compute(),\n",
" small_sky_box_filter.margin.get_partition(order, pixel).compute(),\n",
" ],\n",
" order,\n",
" pixel,\n",
" [\"green\", \"red\"],\n",
" [\"ra\", \"ra\"],\n",
" [\"dec\", \"dec\"],\n",
" xlim=[179.5, 180.1],\n",
" ylim=[9.4, 10.0],\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## Avoiding Duplicates\n",
"\n",
"Joining the margin cache to the catalog's data would introduce duplicate points, where points near the boundary would appear in the margin of one partition, and the main file of another partition. To avoid this, we keep two separate task graphs, one for the catalog and one for its margin. For operations that don't require the margin, the task graphs are kept separate, and when `compute` is called on the catalog, only the catalog's task graph is computed without joining the margin or even loading the margin files. For operations like crossmatching that require the margin, the task graphs are combined with the margin joined and used. For these operations, we use only the margin for the catalog on the right side of the operation. This means that for each left catalog point that is considered, all of the possible nearby matches in the right catalog are also loaded, and so the results are kept accurate. But since there are no duplicates of the left catalog points, there are no duplicate results."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python"
}
Expand Down

0 comments on commit ebcabd1

Please sign in to comment.