examples/GLIMS_ATL06_Subsetter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ca561032-e4a4-4b7a-b01f-46048a33f6e7",
   "metadata": {},
   "source": [
    "## Subset ATL06 to GLIMS Shapefile\n",
    "\n",
    "This notebook uses SlideRule to retrieve ATL06 segments that intersect a provided shapefile.\n",
    "1. Generate the convex hull of the region of interest characterized by the shapefile so that we have a polygon to submit to CMR to get all ATL06 granules that could possibly intersect the region of interest.\n",
    "2. Convert the shapefile to a geojson, and in the process, buffer out the polygons so that no points are missed by SlideRule\n",
    "3. Make the processing request to SlideRule to retrieve all ATL06 segments within the region of interest\n",
    "4. Trim the returned values to the original shapefile to get rid of any segments that were only included in the bufferred region\n",
    "\n",
    "### Notes\n",
    "\n",
    "* SlideRule v4.6.2 and earlier versions have a bug in the ATL06 subsetter that does not handle geojson subsetting correctly.  When running against a SlideRule cluster with this bug, a large amount of data is returned (which takes a long time), and is then substantially trimmed in the final steps.  When running against a later version of SlideRule, the data returned from the server is substantially less and takes significantly less time.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c841bde-7ece-492b-b749-bfc5c9397dbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sliderule import sliderule, icesat2, earthdata\n",
    "from shapely.geometry import Polygon, MultiPolygon, mapping\n",
    "import geopandas as gpd\n",
    "import geojson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a6d95a3-de3b-418a-b2db-62e77887e202",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#\n",
    "# If *.shx file is missing, run this cell to generate\n",
    "#\n",
    "#import fiona\n",
    "#with fiona.Env(SHAPE_RESTORE_SHX='YES'):\n",
    "#    region = sliderule.toregion(\"glims_polygons.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdbb24ce-841c-44aa-a2c7-0c20c2069245",
   "metadata": {},
   "source": [
    "### Read in shapefile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5ac416c-452e-4a04-a578-eae34fda5959",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read shapefile\n",
    "gdf = gpd.read_file(\"glims_polygons.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "114393da-6e0e-4604-9f8e-01c8aa83332d",
   "metadata": {},
   "source": [
    "### Get granules that intersect larger area of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d937c909-3cb3-4e63-b159-2683b3069acd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a multipolygon with simplified internal polygons (needed to get convex hull)\n",
    "polygons = list(gdf.geometry)\n",
    "cleaned_polygons = [polygon.convex_hull for polygon in polygons]\n",
    "cleaned_multipoly = MultiPolygon(cleaned_polygons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe62d6e4-66bd-4ab0-8bd6-d24fcd74deb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# build geojson of multipolygon\n",
    "cleaned_glims_geojson = \"cleaned_glims.geojson\"\n",
    "geojson_obj = geojson.Feature(geometry=mapping(cleaned_multipoly))\n",
    "with open(cleaned_glims_geojson, \"w\") as geojson_file:\n",
    "    geojson.dump(geojson_obj, geojson_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "924fcb89-2c2e-4389-8573-ba528eb6e64d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# get sliderule region of geojson\n",
    "region = sliderule.toregion(cleaned_glims_geojson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87f24364-9e76-4459-a906-9310a5ef0712",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# query CMR for granules that intersect larger area of interest\n",
    "cmr_parms = {\n",
    "    \"asset\": \"icesat2-atl06\",\n",
    "    \"poly\": region[\"poly\"]\n",
    "}\n",
    "earthdata.set_max_resources(350)\n",
    "granules = earthdata.search(cmr_parms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d75d509-7e1b-4730-821a-db855e7df300",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "granules"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "235303b9-ebac-4bf5-90df-c3a2a0011885",
   "metadata": {},
   "source": [
    "### Get detailed geojson of area of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "030ef8a6-8dc8-43cd-8272-66aae250404a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a multipolygon of internal polygons\n",
    "multipoly = MultiPolygon(list(gdf.geometry))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1b2745c-601a-455d-95e8-e4c2df9690ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# buffer out multiplygon\n",
    "buffered_multipoly = multipoly.buffer(0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e68d5f-7cc9-44da-82e1-19317776b7b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# build geojson of multipolygon\n",
    "glims_geojson = \"glims.geojson\"\n",
    "geojson_obj = geojson.Feature(geometry=mapping(buffered_multipoly))\n",
    "with open(glims_geojson, \"w\") as geojson_file:\n",
    "    geojson.dump(geojson_obj, geojson_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "231320d1-2266-4300-a43a-749b4f8a2a89",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = gpd.read_file(\"glims.geojson\")\n",
    "g.plot(markersize=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "684aa4f8-badd-47cc-b7b4-f743c23646d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# open the geojson and read in as raw bytes\n",
    "with open(glims_geojson, mode='rt') as file:\n",
    "    datafile = file.read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1e64edc-57cb-4f15-b812-cca4667b727c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# build the raster parameters for sliderule\n",
    "cellsize = 0.001\n",
    "raster = {\n",
    "    \"data\": datafile, # geojson file\n",
    "    \"length\": len(datafile), # geojson file length\n",
    "    \"cellsize\": cellsize  # units are in crs/projection\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d088011-8858-47d0-89ba-747719397365",
   "metadata": {},
   "source": [
    "### Use sliderule to generate subsetted ATL06 over area of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1518295c-f429-4656-93b8-00f52c52d076",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize the client\n",
    "sliderule.init(\"slideruleearth.io\", verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "659797d8-7991-4684-a64f-d8fe4a12b457",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# atl06 subsetting parameters\n",
    "atl06_parms = {\n",
    "    \"poly\": region[\"poly\"],\n",
    "    \"raster\": raster,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07bcd282-ff63-448e-921d-e924e41343ee",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# make processing request\n",
    "atl06 = icesat2.atl06sp(atl06_parms, resources=granules)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d133f19-5006-4bc8-90a4-fae1f348b54d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# display results\n",
    "atl06"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9c38697-d5a0-4761-b969-fac87d1a996e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot results\n",
    "atl06.plot(markersize=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "992d200b-9951-44d0-8eea-313274fbbb2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save results to a geoparquet file\n",
    "atl06.to_parquet(\"glims_atl06.geoparquet\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d889e500-1309-45ae-b267-de8e4acc78af",
   "metadata": {},
   "source": [
    "### Trim the output to GLIMS polygons\n",
    "The subsetting on SlideRule used a buffered multipolygon so that it wouldn't miss any data.  The steps below further trim the data to the exact region of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c093b5b-9b03-45f5-b6a1-7d7c815d1c26",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read data from geoparquet file, set ICESat-2 crs\n",
    "atl06rb = gpd.read_parquet(\"glims_atl06.geoparquet\")\n",
    "gdf = gdf.set_crs(\"EPSG:7912\", allow_override=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4943490e-69f3-4506-93bc-d31e5357ec91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# trim geodataframe to initial shapefile\n",
    "trimmed_gdf = gpd.sjoin(atl06rb, gdf, how='inner', predicate='within')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea8fcf0e-9f46-422d-b6ce-482b187b2f2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot trimmed results\n",
    "trimmed_gdf.plot(markersize=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "849915e4-927a-41f7-a8b3-9a5de8297bda",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save trimmed results\n",
    "trimmed_gdf.to_parquet(\"glims_subsetted_atl06.geoparquet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77ca35d4-d510-41ed-8c5b-018b77d451ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# display trimmed results\n",
    "trimmed_gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0029f4c-ac1e-4985-b73b-01903132c73c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}