merging latest from development

Author: jpswinski

environment.yml

@@ -6,9 +6,9 @@ dependencies:
   - geopandas
   - h5py
   - ipykernel
-  - ipywidgets>=7.6,<8.0
-  - ipyleaflet>=0.15
-  - jupyterlab=3
+  - ipywidgets
+  - ipyleaflet
+  - jupyterlab
   - jupyterlab_widgets
   - matplotlib
   - ipympl
@@ -25,5 +25,6 @@ dependencies:
   - shapely
   - tk
   - xyzservices
+  - pyarrow
   - pip:
       - -e ./

examples/arcticdem.ipynb

@@ -0,0 +1,148 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "acb12a75-1636-471a-9649-48a408801d4f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from sliderule import icesat2\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib\n",
+    "import geopandas"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b8167cbe-3fe3-4dc9-a5ad-0cbba51c8a07",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "icesat2.init(\"slideruleearth.io\", verbose=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "98ef750a-e88b-4125-b951-d1e29ce50ce2",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "asset = \"nsidc-s3\"\n",
+    "resource = \"ATL03_20190314093716_11600203_005_01.h5\"\n",
+    "region = icesat2.toregion(\"../tests/data/dicksonfjord.geojson\")\n",
+    "parms = { \"poly\": region['poly'],\n",
+    "          \"cnf\": \"atl03_high\",\n",
+    "          \"ats\": 5.0,\n",
+    "          \"cnt\": 5,\n",
+    "          \"len\": 20.0,\n",
+    "          \"res\": 10.0,\n",
+    "          \"maxi\": 1,\n",
+    "          \"samples\": [\"arcticdem-mosaic\"] }\n",
+    "gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d0ebba64-93ab-45c8-9c53-cc5078332617",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gdf[\"delta\"] = gdf[\"h_mean\"] - gdf[\"arcticdem-mosaic-1980-01-06\"]\n",
+    "gdf[\"delta\"].describe()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "db8b0b39-b421-46c6-bb5f-7abc6b1168b7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Setup Plot\n",
+    "fig,ax = plt.subplots(num=None, figsize=(10, 8))\n",
+    "fig.set_facecolor('white')\n",
+    "fig.canvas.header_visible = False\n",
+    "ax.set_title(\"SlideRule vs. ArcticDEM Elevations\")\n",
+    "ax.set_xlabel('UTC')\n",
+    "ax.set_ylabel('height (m)')\n",
+    "legend_elements = []\n",
+    "\n",
+    "# Plot SlideRule ATL06 Elevations\n",
+    "df = gdf[(gdf['rgt'] == 1160) & (gdf['gt'] == 10) & (gdf['cycle'] == 2)]\n",
+    "sc1 = ax.scatter(df.index.values, df[\"h_mean\"].values, c='red', s=2.5)\n",
+    "legend_elements.append(matplotlib.lines.Line2D([0], [0], color='red', lw=6, label='ATL06-SR'))\n",
+    "\n",
+    "# Plot ArcticDEM Elevations\n",
+    "sc2 = ax.scatter(df.index.values, df[\"arcticdem-mosaic-1980-01-06\"].values, c='blue', s=2.5)\n",
+    "legend_elements.append(matplotlib.lines.Line2D([0], [0], color='blue', lw=6, label='ArcticDEM'))\n",
+    "\n",
+    "# Display Legend\n",
+    "lgd = ax.legend(handles=legend_elements, loc=3, frameon=True)\n",
+    "lgd.get_frame().set_alpha(1.0)\n",
+    "lgd.get_frame().set_edgecolor('white')\n",
+    "\n",
+    "# Show Plot\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "82c65e28-468e-463e-9afe-2b52064e7bae",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Setup Plot\n",
+    "fig,ax = plt.subplots(num=None, figsize=(10, 8))\n",
+    "fig.set_facecolor('white')\n",
+    "fig.canvas.header_visible = False\n",
+    "ax.set_title(\"Delta Elevations between SlideRule and ArcticDEM\")\n",
+    "ax.set_xlabel('UTC')\n",
+    "ax.set_ylabel('height (m)')\n",
+    "ax.yaxis.grid(True)\n",
+    "\n",
+    "# Plot Deltas\n",
+    "df = gdf[(gdf['rgt'] == 1160) & (gdf['gt'] == 10) & (gdf['cycle'] == 2)]\n",
+    "sc1 = ax.scatter(df.index.values, df[\"delta\"].values, c='green', s=2.5)\n",
+    "\n",
+    "# Show Plot\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eed8f243-dd0c-4473-a952-fcb2bb863e3c",
+   "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.8.15"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

requirements.txt

@@ -2,5 +2,5 @@ requests
 numpy
 fiona
 geopandas
-shapely<2
+shapely
 scikit-learn

sliderule/icesat2.py

@@ -41,7 +41,6 @@
 import geopandas
 from shapely.geometry.multipolygon import MultiPolygon
 from shapely.geometry import Polygon
-from sklearn.cluster import KMeans
 import sliderule
 
 ###############################################################################
@@ -51,6 +50,14 @@
 # create logger
 logger = logging.getLogger(__name__)
 
+# import cluster support
+clustering_enabled = False
+try:
+    from sklearn.cluster import KMeans
+    clustering_enabled = True
+except:
+    logger.warning("Unable to import sklearn... clustering support disabled")
+
 # profiling times for each major function
 profiles = {}
 
@@ -503,12 +510,8 @@ def __gdf2poly(gdf):
 #
 def __procoutputfile(parm, lon_key, lat_key):
     if "open_on_complete" in parm["output"] and parm["output"]["open_on_complete"]:
-        # Read Parquet File as DataFrame
-        df = geopandas.pd.read_parquet(parm["output"]["path"])
-        # Build GeoDataFrame
-        gdf = __todataframe(df, lon_key=lon_key, lat_key=lat_key)
-        # Return Results
-        return gdf
+        # Return GeoParquet File as GeoDataFrame
+        return geopandas.read_parquet(parm["output"]["path"])
     else:
         # Return Parquet Filename
         return parm["output"]["path"]
@@ -804,7 +807,7 @@ def atl06p(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, callb
             columns = {}
             elevation_records = []
             num_elevations = 0
-            field_dictionary = {} # ['field_name'] = {"extent_id": [], field_name: []}
+            field_dictionary = {} # [<field_name>] = {"extent_id": [], <field_name>: []}
             if len(rsps) > 0:
                 # Sort Records
                 for rsp in rsps:
@@ -814,15 +817,23 @@ def atl06p(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, callb
                     elif 'extrec' == rsp['__rectype']:
                         field_name = parm['atl03_geo_fields'][rsp['field_index']]
                         if field_name not in field_dictionary:
-                            field_dictionary[field_name] = {"extent_id": [], field_name: []}
+                            field_dictionary[field_name] = {'extent_id': [], field_name: []}
                         # Parse Ancillary Data
-                        data = __get_values(rsp['data'], rsp['data_type'], len(rsp['data']))
+                        data = __get_values(rsp['data'], rsp['datatype'], len(rsp['data']))
                         # Add Left Pair Track Entry
                         field_dictionary[field_name]['extent_id'] += rsp['extent_id'] | 0x2,
                         field_dictionary[field_name][field_name] += data[LEFT_PAIR],
                         # Add Right Pair Track Entry
                         field_dictionary[field_name]['extent_id'] += rsp['extent_id'] | 0x3,
                         field_dictionary[field_name][field_name] += data[RIGHT_PAIR],
+                    elif 'rsrec' == rsp['__rectype']:
+                        for sample in rsp["samples"]:
+                            time_str = sliderule.gps2utc(sample["time"])
+                            field_name = parm['samples'][rsp['raster_index']] + "-" + time_str.split(" ")[0].strip()
+                            if field_name not in field_dictionary:
+                                field_dictionary[field_name] = {'extent_id': [], field_name: []}
+                            field_dictionary[field_name]['extent_id'] += rsp['extent_id'],
+                            field_dictionary[field_name][field_name] += sample['value'],
                 # Build Elevation Columns
                 if num_elevations > 0:
                     # Initialize Columns
@@ -967,23 +978,23 @@ def atl03sp(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, call
                         # Get Field Type
                         field_name = parm['atl03_geo_fields'][rsp['field_index']]
                         if field_name not in extent_field_types:
-                            extent_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['data_type']]]["nptype"]
+                            extent_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['datatype']]]["nptype"]
                         # Initialize Extent Dictionary Entry
                         if extent_id not in extent_dictionary:
                             extent_dictionary[extent_id] = {}
                         # Save of Values per Extent ID per Field Name
-                        data = __get_values(rsp['data'], rsp['data_type'], len(rsp['data']))
+                        data = __get_values(rsp['data'], rsp['datatype'], len(rsp['data']))
                         extent_dictionary[extent_id][field_name] = data
                     elif 'phrec' == rsp['__rectype']:
                         # Get Field Type
                         field_name = parm['atl03_ph_fields'][rsp['field_index']]
                         if field_name not in photon_field_types:
-                            photon_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['data_type']]]["nptype"]
+                            photon_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['datatype']]]["nptype"]
                         # Initialize Extent Dictionary Entry
                         if extent_id not in photon_dictionary:
                             photon_dictionary[extent_id] = {}
                         # Save of Values per Extent ID per Field Name
-                        data = __get_values(rsp['data'], rsp['data_type'], len(rsp['data']))
+                        data = __get_values(rsp['data'], rsp['datatype'], len(rsp['data']))
                         photon_dictionary[extent_id][field_name] = data
                 # Build Elevation Columns
                 if num_photons > 0:
@@ -1331,22 +1342,25 @@ def toregion(source, tolerance=0.0, cellsize=0.01, n_clusters=1):
         # generate clusters
         clusters = []
         if n_clusters > 1:
-            # pull out centroids of each geometry object
-            if "CenLon" in gdf and "CenLat" in gdf:
-                X = numpy.column_stack((gdf["CenLon"], gdf["CenLat"]))
+            if clustering_enabled:
+                # pull out centroids of each geometry object
+                if "CenLon" in gdf and "CenLat" in gdf:
+                    X = numpy.column_stack((gdf["CenLon"], gdf["CenLat"]))
+                else:
+                    s = gdf.centroid
+                    X = numpy.column_stack((s.x, s.y))
+                # run k means clustering algorithm against polygons in gdf
+                kmeans = KMeans(n_clusters=n_clusters, init='k-means++', random_state=5, max_iter=400)
+                y_kmeans = kmeans.fit_predict(X)
+                k = geopandas.pd.DataFrame(y_kmeans, columns=['cluster'])
+                gdf = gdf.join(k)
+                # build polygon for each cluster
+                for n in range(n_clusters):
+                    c_gdf = gdf[gdf["cluster"] == n]
+                    c_poly = __gdf2poly(c_gdf)
+                    clusters.append(c_poly)
             else:
-                s = gdf.centroid
-                X = numpy.column_stack((s.x, s.y))
-            # run k means clustering algorithm against polygons in gdf
-            kmeans = KMeans(n_clusters=n_clusters, init='k-means++', random_state=5, max_iter=400)
-            y_kmeans = kmeans.fit_predict(X)
-            k = geopandas.pd.DataFrame(y_kmeans, columns=['cluster'])
-            gdf = gdf.join(k)
-            # build polygon for each cluster
-            for n in range(n_clusters):
-                c_gdf = gdf[gdf["cluster"] == n]
-                c_poly = __gdf2poly(c_gdf)
-                clusters.append(c_poly)
+                raise sliderule.FatalError("Clustering support not enabled; unable to import sklearn package")
 
     # update timing profiles
     profiles[toregion.__name__] = time.perf_counter() - tstart