diff --git a/requirements.txt b/requirements.txt index ce40662..ef505bd 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1,3 @@ xarray==2023.8.0 -odc-stac \ No newline at end of file +odc-stac +rasterstats diff --git a/stac-R.ipynb b/stac-R.ipynb index 39daf7d..3bfc9ff 100644 --- a/stac-R.ipynb +++ b/stac-R.ipynb @@ -32,7 +32,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "id": "f8567e1c-e63c-4ce5-8175-5f14b0ae3312", "metadata": { "tags": [], @@ -66,19 +66,21 @@ }, "outputs": [], "source": [ - "col <-\n", - " stac_image_collection(items$features,\n", - " asset_names = c(\"B04\",\"B08\", \"SCL\"),\n", - " property_filter = \\(x) {x[[\"eo:cloud_cover\"]] < 20})\n", "\n", - "cube <- cube_view(srs = \"EPSG:4326\", \n", + "col <- stac_image_collection(items$features,\n", + " asset_names = c(\"B04\", \"B08\", \"SCL\"),\n", + " property_filter = \\(x){\n", + " x[[\"eo:cloud_cover\"]] < 20\n", + " })\n", + "\n", + "cube <- cube_view(srs = \"EPSG:4326\",\n", " extent = list(t0 = start_date, t1 = end_date,\n", " left = box[1], right = box[3],\n", " top = box[4], bottom = box[2]),\n", " nx = 1200, ny = 1200, dt = \"P1D\",\n", " aggregation = \"median\", resampling = \"average\")\n", "\n", - "mask <- image_mask(\"SCL\", values=c(3,8,9)) # mask clouds and cloud shadows\n" + "mask <- image_mask(\"SCL\", values=c(3, 8, 9)) # mask clouds and cloud shadows\n" ] }, { @@ -113,13 +115,28 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": { "vscode": { "languageId": "r" } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading layer `CASanFrancisco1937' from data source \n", + " `/vsicurl/https://dsl.richmond.edu/panorama/redlining/static/downloads/geojson/CASanFrancisco1937.geojson' \n", + " using driver `GeoJSON'\n", + "Simple feature collection with 97 features and 4 fields\n", + "Geometry type: MULTIPOLYGON\n", + "Dimension: XY\n", + "Bounding box: xmin: -122.5101 ymin: 37.70801 xmax: -122.3627 ymax: 37.80668\n", + "Geodetic CRS: WGS 84\n" + ] + } + ], "source": [ "sf <- st_read(\"/vsicurl/https://dsl.richmond.edu/panorama/redlining/static/downloads/geojson/CASanFrancisco1937.geojson\") |>\n", " st_make_valid()\n", @@ -143,13 +160,7 @@ "output_type": "stream", "text": [ "Warning message:\n", - "“v3 code detected: as of tmap v4, tm_legend should be specified per visual variable (e.g. with the argument fill.legend of tm_polygons”\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ + "“v3 code detected: as of tmap v4, tm_legend should be specified per visual variable (e.g. with the argument fill.legend of tm_polygons”\n", "Warning message in value[[3L]](cond):\n", "“could not rename the data.table”\n" ] diff --git a/stac-odc.ipynb b/stac-odc.ipynb index 39776c7..ce3b6aa 100644 --- a/stac-odc.ipynb +++ b/stac-odc.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, + "execution_count": 14, "id": "ae1766b4-e3b0-414c-a12e-d07f4d001a7e", "metadata": { "tags": [] @@ -20,7 +20,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 15, "id": "2fa33351-ed8f-4698-8bfc-2fbdb413d7cc", "metadata": { "tags": [] @@ -39,7 +39,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 16, "id": "db5e1558-d227-4a13-8ed7-1fed6e21e225", "metadata": { "tags": [] @@ -52,7 +52,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 17, "metadata": {}, "outputs": [], "source": [ @@ -68,23 +68,16 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 18, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "computing ndvi...\n" - ] - } - ], + "outputs": [], "source": [ "\n", "red = data.red\n", "nir = data.nir08\n", "\n", - "print(\"computing ndvi...\")\n", + "# summarize over time. \n", + "# quite probably better to use resampling strategy in odc.stac.load though.\n", "import dask.diagnostics\n", "with dask.diagnostics.ProgressBar():\n", " ndvi = ( ((nir - red) / (red + nir)).\n", @@ -94,21 +87,21 @@ " )\n", "\n", "# mask out bad pixels\n", - "mask = ndvi.where(ndvi <= 1)\n" + "ndvi = ndvi.where(ndvi <= 1)" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 10, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" }, @@ -129,8 +122,70 @@ "cmap = plt.colormaps.get_cmap('viridis') # viridis is the default colormap for imshow\n", "cmap.set_bad(color='black')\n", "\n", + "ndvi.plot.imshow(row=\"time\", cmap=cmap, add_colorbar=False, size=4)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import geopandas as gpd\n", + "sf = \"/vsicurl/https://dsl.richmond.edu/panorama/redlining/static/downloads/geojson/CASanFrancisco1937.geojson\"\n", + "poly = gpd.read_file(sf)\n", + "poly = poly.to_crs(epsg=\"32610\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "xarray.core.dataarray.DataArray" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "Specify affine transform for numpy arrays", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m/workspaces/nasa-topst-env-justice/stac-odc.ipynb Cell 9\u001b[0m line \u001b[0;36m3\n\u001b[1;32m 1\u001b[0m \u001b[39mimport\u001b[39;00m \u001b[39mrasterstats\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m result \u001b[39m=\u001b[39m rasterstats\u001b[39m.\u001b[39;49mzonal_stats(\n\u001b[1;32m 4\u001b[0m vectors \u001b[39m=\u001b[39;49m sf, \n\u001b[1;32m 5\u001b[0m raster \u001b[39m=\u001b[39;49m ndvi\u001b[39m.\u001b[39;49mvalues,\n\u001b[1;32m 6\u001b[0m stats \u001b[39m=\u001b[39;49m [\u001b[39m'\u001b[39;49m\u001b[39mmean\u001b[39;49m\u001b[39m'\u001b[39;49m, \u001b[39m'\u001b[39;49m\u001b[39mmin\u001b[39;49m\u001b[39m'\u001b[39;49m, \u001b[39m'\u001b[39;49m\u001b[39mmax\u001b[39;49m\u001b[39m'\u001b[39;49m]\n\u001b[1;32m 7\u001b[0m )\n\u001b[1;32m 8\u001b[0m result\n", + "File \u001b[0;32m~/.local/lib/python3.10/site-packages/rasterstats/main.py:36\u001b[0m, in \u001b[0;36mzonal_stats\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 28\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mzonal_stats\u001b[39m(\u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs):\n\u001b[1;32m 29\u001b[0m \u001b[39m \u001b[39m\u001b[39m\"\"\"The primary zonal statistics entry point.\u001b[39;00m\n\u001b[1;32m 30\u001b[0m \n\u001b[1;32m 31\u001b[0m \u001b[39m All arguments are passed directly to ``gen_zonal_stats``.\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 34\u001b[0m \u001b[39m The only difference is that ``zonal_stats`` will\u001b[39;00m\n\u001b[1;32m 35\u001b[0m \u001b[39m return a list rather than a generator.\"\"\"\u001b[39;00m\n\u001b[0;32m---> 36\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mlist\u001b[39;49m(gen_zonal_stats(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs))\n", + "File \u001b[0;32m~/.local/lib/python3.10/site-packages/rasterstats/main.py:161\u001b[0m, in \u001b[0;36mgen_zonal_stats\u001b[0;34m(vectors, raster, layer, band, nodata, affine, stats, all_touched, categorical, category_map, add_stats, zone_func, raster_out, prefix, geojson_out, boundless, **kwargs)\u001b[0m\n\u001b[1;32m 158\u001b[0m warnings\u001b[39m.\u001b[39mwarn(\u001b[39m\"\u001b[39m\u001b[39mUse `band` to specify band number\u001b[39m\u001b[39m\"\u001b[39m, \u001b[39mDeprecationWarning\u001b[39;00m)\n\u001b[1;32m 159\u001b[0m band \u001b[39m=\u001b[39m band_num\n\u001b[0;32m--> 161\u001b[0m \u001b[39mwith\u001b[39;00m Raster(raster, affine, nodata, band) \u001b[39mas\u001b[39;00m rast:\n\u001b[1;32m 162\u001b[0m features_iter \u001b[39m=\u001b[39m read_features(vectors, layer)\n\u001b[1;32m 163\u001b[0m \u001b[39mfor\u001b[39;00m _, feat \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(features_iter):\n", + "File \u001b[0;32m~/.local/lib/python3.10/site-packages/rasterstats/io.py:261\u001b[0m, in \u001b[0;36mRaster.__init__\u001b[0;34m(self, raster, affine, nodata, band)\u001b[0m\n\u001b[1;32m 259\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(raster, np\u001b[39m.\u001b[39mndarray):\n\u001b[1;32m 260\u001b[0m \u001b[39mif\u001b[39;00m affine \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m--> 261\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mSpecify affine transform for numpy arrays\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m 262\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39marray \u001b[39m=\u001b[39m raster\n\u001b[1;32m 263\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39maffine \u001b[39m=\u001b[39m affine\n", + "\u001b[0;31mValueError\u001b[0m: Specify affine transform for numpy arrays" + ] + } + ], + "source": [ + "import rasterstats\n", "\n", - "mask.plot.imshow(row=\"time\", cmap=cmap, add_colorbar=False, size=4)\n" + "result = rasterstats.zonal_stats(\n", + " vectors = sf, \n", + " raster = ndvi.values,\n", + " nodata = src_srtm.nodata, \n", + " affine = src_srtm.transform,\n", + " stats = ['mean', 'min', 'max']\n", + ")\n", + "result" ] } ], diff --git a/stac.ipynb b/stac.ipynb index 1b00c4f..4a45b55 100644 --- a/stac.ipynb +++ b/stac.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "id": "ae1766b4-e3b0-414c-a12e-d07f4d001a7e", "metadata": { "tags": [] @@ -18,7 +18,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "id": "2fa33351-ed8f-4698-8bfc-2fbdb413d7cc", "metadata": { "tags": [] @@ -37,7 +37,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "id": "db5e1558-d227-4a13-8ed7-1fed6e21e225", "metadata": { "tags": [] @@ -409,7 +409,7 @@ " stroke: currentColor;\n", " fill: currentColor;\n", "}\n", - "
<xarray.DataArray 'stackstac-32caef4f1496167f7ad4b58bb5aacca8' (time: 3,\n",
+       "
<xarray.DataArray 'stackstac-85cc2488c6f5b7f5a7bb46662333113c' (time: 3,\n",
        "                                                                band: 32,\n",
        "                                                                y: 1103, x: 1306)>\n",
        "dask.array<fetch_raster_window, shape=(3, 32, 1103, 1306), dtype=float64, chunksize=(1, 1, 1024, 1024), chunktype=numpy.ndarray>\n",
@@ -419,7 +419,7 @@
        "  * band                                     (band) <U12 'aot' ... 'wvp-jp2'\n",
        "  * x                                        (x) float64 5.431e+05 ... 5.562e+05\n",
        "  * y                                        (y) float64 4.185e+06 ... 4.174e+06\n",
-       "    s2:granule_id                            (time) <U62 'S2A_OPER_MSI_L2A_TL...\n",
+       "    s2:not_vegetated_percentage              (time) float64 23.64 23.44 24.6\n",
        "    ...                                       ...\n",
        "    raster:bands                             (band) object [{'nodata': 0, 'da...\n",
        "    gsd                                      (band) object None 10 ... None None\n",
@@ -431,7 +431,7 @@
        "    spec:        RasterSpec(epsg=32610, bounds=(543120.0, 4173530.0, 556180.0...\n",
        "    crs:         epsg:32610\n",
        "    transform:   | 10.00, 0.00, 543120.00|\\n| 0.00,-10.00, 4184560.00|\\n| 0.0...\n",
-       "    resolution:  10.0
  • spec :
    RasterSpec(epsg=32610, bounds=(543120.0, 4173530.0, 556180.0, 4184560.0), resolutions_xy=(10.0, 10.0))
    crs :
    epsg:32610
    transform :
    | 10.00, 0.00, 543120.00|\n", "| 0.00,-10.00, 4184560.00|\n", "| 0.00, 0.00, 1.00|
    resolution :
    10.0
    • " ], "text/plain": [ - "\n", "dask.array\n", @@ -712,7 +712,7 @@ " * band (band) " + "" ] }, - "execution_count": 25, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" },