updated raster support for zonal statistics

Author: jpswinski

examples/arcticdem.ipynb

@@ -42,25 +42,57 @@
     "          \"len\": 20.0,\n",
     "          \"res\": 10.0,\n",
     "          \"maxi\": 1,\n",
-    "          \"samples\": [\"arcticdem-mosaic\"] }\n",
+    "          \"samples\": {\"mosaic\": {\"asset\": \"arcticdem-mosaic\", \"radius\": 10.0, \"zonal_stats\": True}} }\n",
     "gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "d0ebba64-93ab-45c8-9c53-cc5078332617",
+   "id": "fd15bf14-ab10-4cf6-9524-592962a8f8b2",
    "metadata": {},
    "outputs": [],
    "source": [
-    "gdf[\"delta\"] = gdf[\"h_mean\"] - gdf[\"arcticdem-mosaic-1980-01-06\"]\n",
-    "gdf[\"delta\"].describe()"
+    "gdf"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "db8b0b39-b421-46c6-bb5f-7abc6b1168b7",
+   "id": "59ea096b-3443-4b9a-b114-e818f143ca45",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gdf[\"value_delta\"] = gdf[\"h_mean\"] - gdf[\"mosaic.value\"]\n",
+    "gdf[\"value_delta\"].describe()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f80cd750-9b91-406a-ac6a-dc04937cd811",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gdf[\"mean_delta\"] = gdf[\"h_mean\"] - gdf[\"mosaic.mean\"]\n",
+    "gdf[\"mean_delta\"].describe()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4e8fe135-bc15-4b8b-887f-ae16ac81487f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gdf[\"median_delta\"] = gdf[\"h_mean\"] - gdf[\"mosaic.median\"]\n",
+    "gdf[\"median_delta\"].describe()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d6785ed8-cb0d-49cb-9de1-30ac5792884a",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -79,9 +111,13 @@
     "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",
+    "sc2 = ax.scatter(df.index.values, df[\"mosaic.value\"].values, c='blue', s=2.5)\n",
     "legend_elements.append(matplotlib.lines.Line2D([0], [0], color='blue', lw=6, label='ArcticDEM'))\n",
     "\n",
+    "# Plot ArcticDEM Mean Elevations\n",
+    "sc2 = ax.scatter(df.index.values, df[\"mosaic.value\"].values, c='green', s=2.5)\n",
+    "legend_elements.append(matplotlib.lines.Line2D([0], [0], color='green', 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",
@@ -108,8 +144,12 @@
     "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",
+    "df1 = gdf[(gdf['rgt'] == 1160) & (gdf['gt'] == 10) & (gdf['cycle'] == 2)]\n",
+    "sc1 = ax.scatter(df1.index.values, df1[\"value_delta\"].values, c='blue', s=2.5)\n",
+    "\n",
+    "# Plot Deltas\n",
+    "df2 = gdf[(gdf['rgt'] == 1160) & (gdf['gt'] == 10) & (gdf['cycle'] == 2)]\n",
+    "sc2 = ax.scatter(df2.index.values, df2[\"mean_delta\"].values, c='green', s=2.5)\n",
     "\n",
     "# Show Plot\n",
     "plt.show()"

sliderule/icesat2.py

@@ -826,14 +826,33 @@ def atl06p(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, callb
                         # 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'],
+                    elif 'rsrec' == rsp['__rectype'] or 'zsrec' == rsp['__rectype']:
+                        if rsp["num_samples"] <= 0:
+                            continue
+                        # Get field names and set
+                        sample = rsp["samples"][0]
+                        field_names = list(sample.keys())
+                        field_names.remove("__rectype")
+                        field_set = rsp['key']
+                        as_numpy_array = False
+                        if rsp["num_samples"] > 1:
+                            as_numpy_array = True
+                        # On first time, build empty dictionary for field set associated with raster
+                        if field_set not in field_dictionary:
+                            field_dictionary[field_set] = {'extent_id': []}
+                            for field in field_names:
+                                field_dictionary[field_set][field_set + "." + field] = []
+                        # Populate dictionary for field set
+                        field_dictionary[field_set]['extent_id'] += rsp['extent_id'],
+                        for field in field_names:
+                            if as_numpy_array:
+                                data = []
+                                for s in rsp["samples"]:
+                                    data += s[field],
+                                field_dictionary[field_set][field_set + "." + field] += numpy.array(data),
+                            else:
+                                field_dictionary[field_set][field_set + "." + field] += sample[field],
+
                 # Build Elevation Columns
                 if num_elevations > 0:
                     # Initialize Columns

tests/test_arcticdem.py

@@ -17,7 +17,7 @@ def test_vrt(self, domain, organization):
         assert abs(rsps["samples"][0][0]["value"] - 80.713500976562) < 0.001
         assert rsps["samples"][0][0]["file"] == '/vsis3/pgc-opendata-dems/arcticdem/mosaics/v3.0/2m/70_09/70_09_2_1_2m_v3.0_reg_dem.tif'
 
-    def test_sampler(self, domain, asset, organization):
+    def test_nearestneighbour(self, domain, asset, organization):
         icesat2.init(domain, organization=organization)
         resource = "ATL03_20190314093716_11600203_005_01.h5"
         region = icesat2.toregion(os.path.join(TESTDIR, "data/dicksonfjord.geojson"))
@@ -29,12 +29,37 @@ def test_sampler(self, domain, asset, organization):
                   "len": 40.0,
                   "res": 20.0,
                   "maxi": 1,
-                  "samples": ["arcticdem-mosaic"] }
+                  "samples": {"mosaic": {"asset": "arcticdem-mosaic"}} }
         gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])
         assert len(gdf) == 964
-        assert len(gdf.keys()) == 17
+        assert len(gdf.keys()) == 18
         assert gdf["rgt"][0] == 1160
         assert gdf["cycle"][0] == 2
         assert gdf['segment_id'].describe()["min"] == 405240
         assert gdf['segment_id'].describe()["max"] == 405915
-        assert abs(gdf["arcticdem-mosaic-1980-01-06"].describe()["min"] - 655.14990234375) < 0.0001
+        assert abs(gdf["mosaic.value"].describe()["min"] - 655.14990234375) < 0.0001
+
+    def test_zonal_stats(self, domain, asset, organization):
+        icesat2.init(domain, organization=organization)
+        resource = "ATL03_20190314093716_11600203_005_01.h5"
+        region = icesat2.toregion(os.path.join(TESTDIR, "data/dicksonfjord.geojson"))
+        parms = { "poly": region['poly'],
+                  "raster": region['raster'],
+                  "cnf": "atl03_high",
+                  "ats": 20.0,
+                  "cnt": 10,
+                  "len": 40.0,
+                  "res": 20.0,
+                  "maxi": 1,
+                  "samples": {"mosaic": {"asset": "arcticdem-mosaic", "radius": 10.0, "zonal_stats": True}} }
+        gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])
+        assert len(gdf) == 964
+        assert len(gdf.keys()) == 25
+        assert gdf["rgt"][0] == 1160
+        assert gdf["cycle"][0] == 2
+        assert gdf['segment_id'].describe()["min"] == 405240
+        assert gdf['segment_id'].describe()["max"] == 405915
+        assert abs(gdf["mosaic.value"].describe()["min"] - 655.14990234375) < 0.0001
+        assert gdf["mosaic.count"].describe()["max"] == 81
+        assert gdf["mosaic.stdev"].describe()["count"] == 964
+        assert gdf["mosaic.time"][0] == 1176076818.0