Skip to content

Commit

Permalink
🖊️
Browse files Browse the repository at this point in the history
  • Loading branch information
@cboettig
cboettig committed Nov 16, 2023
1 parent e21e55f commit b015694
Show file tree
Hide file tree
Showing 4 changed files with 146 additions and 79 deletions.
3 changes: 2 additions & 1 deletion requirements.txt
View file
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
xarray==2023.8.0
odc-stac
odc-stac
rasterstats
43 changes: 27 additions & 16 deletions stac-R.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 2,
"id": "f8567e1c-e63c-4ce5-8175-5f14b0ae3312",
"metadata": {
"tags": [],
Expand Down Expand Up @@ -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"
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -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"
]
Expand Down
95 changes: 75 additions & 20 deletions stac-odc.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 14,
"id": "ae1766b4-e3b0-414c-a12e-d07f4d001a7e",
"metadata": {
"tags": []
Expand All @@ -20,7 +20,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 15,
"id": "2fa33351-ed8f-4698-8bfc-2fbdb413d7cc",
"metadata": {
"tags": []
Expand All @@ -39,7 +39,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 16,
"id": "db5e1558-d227-4a13-8ed7-1fed6e21e225",
"metadata": {
"tags": []
Expand All @@ -52,7 +52,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -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",
Expand All @@ -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": [
"<xarray.plot.facetgrid.FacetGrid at 0x7fdb22a1a500>"
"<xarray.plot.facetgrid.FacetGrid at 0x7f46d3f07820>"
]
},
"execution_count": 10,
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
},
Expand All @@ -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 <a href='vscode-notebook-cell://codespaces%2Bobscure-spoon-jqgxqq6v3qg96/workspaces/nasa-topst-env-justice/stac-odc.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=0'>1</a>\u001b[0m \u001b[39mimport\u001b[39;00m \u001b[39mrasterstats\u001b[39;00m\n\u001b[0;32m----> <a href='vscode-notebook-cell://codespaces%2Bobscure-spoon-jqgxqq6v3qg96/workspaces/nasa-topst-env-justice/stac-odc.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=2'>3</a>\u001b[0m result \u001b[39m=\u001b[39m rasterstats\u001b[39m.\u001b[39;49mzonal_stats(\n\u001b[1;32m <a href='vscode-notebook-cell://codespaces%2Bobscure-spoon-jqgxqq6v3qg96/workspaces/nasa-topst-env-justice/stac-odc.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=3'>4</a>\u001b[0m vectors \u001b[39m=\u001b[39;49m sf, \n\u001b[1;32m <a href='vscode-notebook-cell://codespaces%2Bobscure-spoon-jqgxqq6v3qg96/workspaces/nasa-topst-env-justice/stac-odc.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=4'>5</a>\u001b[0m raster \u001b[39m=\u001b[39;49m ndvi\u001b[39m.\u001b[39;49mvalues,\n\u001b[1;32m <a href='vscode-notebook-cell://codespaces%2Bobscure-spoon-jqgxqq6v3qg96/workspaces/nasa-topst-env-justice/stac-odc.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=5'>6</a>\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 <a href='vscode-notebook-cell://codespaces%2Bobscure-spoon-jqgxqq6v3qg96/workspaces/nasa-topst-env-justice/stac-odc.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=6'>7</a>\u001b[0m )\n\u001b[1;32m <a href='vscode-notebook-cell://codespaces%2Bobscure-spoon-jqgxqq6v3qg96/workspaces/nasa-topst-env-justice/stac-odc.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=7'>8</a>\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"
]
}
],
Expand Down
Loading

0 comments on commit b015694

Please sign in to comment.