add s3 retrieve and granule plot

Author: tsutterley

examples/cmr_debug_regions.ipynb

@@ -20,18 +20,52 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from sliderule import icesat2, ipysliderule\n",
-    "import numpy as np\n",
+    "import re\n",
+    "import time\n",
     "import logging\n",
+    "import posixpath\n",
+    "import numpy as np\n",
+    "import geopandas as gpd\n",
+    "import concurrent.futures\n",
     "import matplotlib.cm as cm\n",
     "import matplotlib.colors as colors\n",
+    "import ipywidgets as widgets\n",
+    "from shapely.geometry import LineString\n",
+    "from sliderule import sliderule, icesat2, ipysliderule\n",
+    "import sliderule.io\n",
     "# autoreload\n",
     "%load_ext autoreload\n",
     "%autoreload 2\n",
     "# create logger\n",
     "logging.basicConfig(level=logging.INFO)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Set ICESat-2 Product"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Configure ICESat-2 API\n",
+    "icesat2.init(\"icesat2sliderule.org\", loglevel=logging.WARNING)\n",
+    "icesat2.get_version()\n",
+    "\n",
+    "# dropdown menu for ICESat-2 product\n",
+    "product_select = widgets.Dropdown(\n",
+    "    options=['ATL03','ATL06','ATL08'],\n",
+    "    description='Product:',\n",
+    "    disabled=False\n",
+    ")\n",
+    "display(product_select)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -68,7 +102,9 @@
     "granule_polygons = []\n",
     "for poly in m.regions:\n",
     "    # polygon from map\n",
-    "    resources,polygons = icesat2.cmr(polygon=poly, return_polygons=True)\n",
+    "    resources,polygons = icesat2.cmr(polygon=poly,\n",
+    "        short_name=product_select.value,\n",
+    "        return_polygons=True)\n",
     "    granule_list.extend(resources)\n",
     "    granule_polygons.extend(polygons)\n",
     "# print list of granules\n",
@@ -77,6 +113,21 @@
     "logging.debug(granule_list)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "granule_select = widgets.SelectMultiple(\n",
+    "    options=granule_list,\n",
+    "    description='Granules:',\n",
+    "    layout=widgets.Layout(width='35%', height='200px'),\n",
+    "    disabled=False\n",
+    ")\n",
+    "display(granule_select)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -90,9 +141,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cmap = iter(cm.viridis(np.linspace(0,1,num_granules)))\n",
-    "for region in granule_polygons:\n",
-    "    locations = [(p['lat'],p['lon']) for p in region]\n",
+    "granule_indices = list(granule_select.index)\n",
+    "cmap = iter(cm.viridis(np.linspace(0,1,len(granule_indices))))\n",
+    "for g in granule_indices:\n",
+    "    locations = [(p['lat'],p['lon']) for p in granule_polygons[g]]\n",
     "    color = colors.to_hex(next(cmap))\n",
     "    polygon = ipysliderule.ipyleaflet.Polygon(\n",
     "        locations=locations,\n",
@@ -103,6 +155,146 @@
     "    )\n",
     "    m.map.add_layer(polygon)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def s3_retrieve(granule, **kwargs):\n",
+    "    # set default keyword arguments\n",
+    "    kwargs.setdefault('asset','nsidc-s3')\n",
+    "    kwargs.setdefault('lon_key','reference_photon_lon')\n",
+    "    kwargs.setdefault('lat_key','reference_photon_lat')\n",
+    "    kwargs.setdefault('index_key','time')\n",
+    "    kwargs.setdefault('polygon',None)\n",
+    "    # regular expression operator for extracting information from files\n",
+    "    rx = re.compile(r'(ATL\\d{2})(-\\d{2})?_(\\d{4})(\\d{2})(\\d{2})(\\d{2})'\n",
+    "        r'(\\d{2})(\\d{2})_(\\d{4})(\\d{2})(\\d{2})_(\\d{3})_(\\d{2})(.*?).h5$')\n",
+    "    # extract parameters from ICESat-2 granule\n",
+    "    PRD,HEM,YY,MM,DD,HH,MN,SS,TRK,CYCL,GRN,RL,VRS,AUX=rx.findall(granule).pop()\n",
+    "    # variables of interest\n",
+    "    if (PRD == 'ATL03'):\n",
+    "        segment_group = \"geolocation\"\n",
+    "        segment_key = 'segment_id'\n",
+    "        vnames = ['segment_id','delta_time','reference_photon_lat',\n",
+    "            'reference_photon_lon']\n",
+    "    elif (PRD == 'ATL06'):\n",
+    "        segment_group = \"land_ice_segments\"\n",
+    "        segment_key = 'segment_id'\n",
+    "        vnames = ['segment_id','delta_time','latitude','longitude']\n",
+    "    elif (PRD == 'ATL08'):\n",
+    "        segment_group = \"land_segments\"\n",
+    "        segment_key = 'segment_id_beg'\n",
+    "        vnames = ['segment_id_beg','segment_id_end','delta_time',\n",
+    "            'latitude','longitude']\n",
+    "    # for each valid beam within the HDF5 file\n",
+    "    frames = []\n",
+    "    gt = dict(gt1l=10,gt1r=20,gt2l=30,gt2r=40,gt3l=50,gt3r=60)\n",
+    "    atlas_sdp_epoch = np.datetime64('2018-01-01T00:00:00')\n",
+    "    kwds = dict(startrow=0,numrows=-1)\n",
+    "    for gtx in ['gt1l','gt1r','gt2l','gt2r','gt3l','gt3r']:\n",
+    "        geodatasets = [dict(dataset=f'{gtx}/{segment_group}/{v}',**kwds) for v in vnames]\n",
+    "        try:\n",
+    "            # get datasets from s3\n",
+    "            hidatasets = icesat2.h5p(geodatasets, granule, kwargs['asset'])\n",
+    "            # copy to new \"flattened\" dictionary\n",
+    "            data = {posixpath.basename(key):var for key,var in hidatasets.items()}\n",
+    "            # Generate Time Column\n",
+    "            delta_time = (data['delta_time']*1e9).astype('timedelta64[ns]')\n",
+    "            data['time'] = gpd.pd.to_datetime(atlas_sdp_epoch + delta_time)\n",
+    "        except:\n",
+    "            pass\n",
+    "        else:\n",
+    "            # copy filename parameters\n",
+    "            data['rgt'] = [int(TRK)]*len(data['delta_time'])\n",
+    "            data['cycle'] = [int(CYCL)]*len(data['delta_time'])\n",
+    "            data['gt'] = [gt[gtx]]*len(data['delta_time'])\n",
+    "            # pandas dataframe from compiled dictionary\n",
+    "            frames.append(gpd.pd.DataFrame.from_dict(data))\n",
+    "    # concatenate pandas dataframe\n",
+    "    try:\n",
+    "        df = gpd.pd.concat(frames)\n",
+    "    except:\n",
+    "        return sliderule.icesat2.__emptyframe()\n",
+    "    # convert to a GeoDataFrame\n",
+    "    lon_key,lat_key = (kwargs['lon_key'],kwargs['lat_key'])\n",
+    "    geometry = gpd.points_from_xy(df[lon_key], df[lat_key])\n",
+    "    gdf = gpd.GeoDataFrame(df.drop(columns=[lon_key,lat_key]),\n",
+    "        geometry=geometry,crs='EPSG:4326')\n",
+    "    # create global track variable\n",
+    "    track = 6*(gdf['rgt']-1) + (gdf['gt']/10)\n",
+    "    gdf = gdf.assign(track=track)\n",
+    "    # calculate global reference point\n",
+    "    global_ref_pt = 6*1387*gdf[segment_key] + track\n",
+    "    gdf = gdf.assign(global_ref_pt=global_ref_pt)\n",
+    "    # sort values for reproducible output despite async processing\n",
+    "    gdf.set_index(kwargs['index_key'], inplace=True)\n",
+    "    gdf.sort_index(inplace=True)\n",
+    "    # remove duplicate points\n",
+    "    gdf = gdf[~gdf.index.duplicated()]\n",
+    "    # intersect with geometry in projected reference system\n",
+    "    if kwargs['polygon'] is not None:\n",
+    "        gdf = gpd.overlay(gdf.to_crs(kwargs['polygon'].crs),\n",
+    "            kwargs['polygon'], how='intersection')\n",
+    "    # convert back to original coordinate reference system\n",
+    "    return gdf.to_crs('EPSG:4326')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# granule resources for selected segments\n",
+    "perf_start = time.perf_counter()\n",
+    "gdf = sliderule.icesat2.__emptyframe()\n",
+    "num_servers = sliderule.update_available_servers()\n",
+    "max_workers = num_servers * icesat2.SERVER_SCALE_FACTOR\n",
+    "with concurrent.futures.ThreadPoolExecutor(max_workers=max_workers) as executor:\n",
+    "    futures = [executor.submit(s3_retrieve, granule_list[g]) for g in granule_indices]\n",
+    "    # Wait for Results\n",
+    "    for future in concurrent.futures.as_completed(futures):\n",
+    "        # append to dataframe\n",
+    "        gdf = gdf.append(future.result())\n",
+    "\n",
+    "# Display Statistics\n",
+    "perf_stop = time.perf_counter()\n",
+    "perf_duration = perf_stop - perf_start\n",
+    "print(\"Completed in {:.3f} seconds of wall-clock time\".format(perf_duration))\n",
+    "print(\"Reference Ground Tracks: {}\".format(gdf[\"rgt\"].unique()))\n",
+    "print(\"Cycles: {}\".format(gdf[\"cycle\"].unique()))\n",
+    "print(\"Received {} elevations\".format(gdf.shape[0]))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# fix int columns that were converted in objects\n",
+    "fixed = gdf.drop(columns=['geometry'])\n",
+    "for column in fixed.select_dtypes(include='object').columns:\n",
+    "    fixed[column] = fixed[column].astype(\"int32\")\n",
+    "fixed = gpd.GeoDataFrame(fixed,geometry=gdf.geometry,crs='EPSG:4326')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# convert from points to linestrings grouping by track\n",
+    "grouped = fixed.groupby(['track'])['geometry'].apply(\n",
+    "    lambda x: LineString(x.tolist()) if x.size > 1 else x.tolist())\n",
+    "geodata = ipysliderule.ipyleaflet.GeoData(geo_dataframe=gpd.GeoDataFrame(grouped),\n",
+    "    style={'color': 'black', 'opacity':1, 'weight':0.1})\n",
+    "m.map.add_layer(geodata)"
+   ]
   }
  ],
  "metadata": {