Modernize notebook: heat_and_trees

heat_and_trees/anaconda-project.yml

@@ -15,39 +15,28 @@ examples_config:
 
 user_fields: [examples_config]
 
-channels: [defaults]
+# channels: [defaults]
+channels: [conda-forge, nodefaults]
 
 packages: &pkgs
-- python=3.8
-- notebook
-- bokeh <3
-- cartopy <0.21
-- colorcet
-- datashader <0.14
-- fastparquet =0.5.0
-- geopandas =0.9.0
-- geoviews =1.9.1
-- holoviews =1.14.4
-- hvplot =0.7.3
-- intake =0.6.2
-- intake-xarray =0.5.0
-- numpy <1.24
-- pandas <1.4
-- pip
-- rasterio <1.3
-- tqdm =4.61.2
-- xarray =0.19.0
-- param =1.11.1
-- aiohttp =3.7.4
-- requests =2.25.1
-- pyepsg =0.4.0
-- numba <0.56
-- panel =0.12.0
-- jinja2 <3
-- markupSafe =2.0.1
-- shapely <2
-- pip:
-  - rio-toa
+- python=3.11
+- notebook<7
+- bokeh>=3.4.1
+- cartopy>=0.23.0
+- colorcet>=3.1.0
+- datashader>=0.16.1
+- geopandas>=0.14.4
+- geoviews>=1.12.0
+- hvplot>=0.10.0
+- intake<2
+- aiohttp>=3.9.5
+- intake-xarray>=0.7.0
+- numpy>=1.26.4
+- pandas>=2.2.2
+- pyarrow>=16.1.0
+- rasterio>=1.2.10
+- rioxarray>=0.15.0
+- xarray>=2024.5.0
 
 dependencies: *pkgs
 
@@ -58,8 +47,6 @@ commands:
 variables:
   INTAKE_CACHE_DIR: data
 
-downloads: {}
-
 platforms:
 - linux-64
 - osx-64

heat_and_trees/catalog.yml

@@ -24,3 +24,13 @@ sources:
         band: 1
         x: 256
         y: 256
+  neighborhoods_philadelphia:
+    description: Neighborhoods Philadelphia
+    driver: json
+    args:
+      urlpath: 'https://datasets.holoviz.org/philadelphia_trees/v1/Neighborhoods_Philadelphia.geojson'
+  ppr_tree_inventory_2022:
+    description: PPR Tree Inventory 2022
+    driver: json
+    args:
+      urlpath: 'https://datasets.holoviz.org/philadelphia_trees/v1/PPR_Tree_Inventory_2022.geojson'

heat_and_trees/heat_and_trees.ipynb

@@ -12,12 +12,9 @@
    "metadata": {},
    "source": [
     "## Urban Heat Islands and Street Trees\n",
-    "In this notebook we'll be exploring the urban heat island effect by looking at the impact on surface temperature of roof color and street trees. We'll be replicating the process described here: http://urbanspatialanalysis.com/urban-heat-islands-street-trees-in-philadelphia/ but using Python tools rather than ESRI.\n",
+    "In this notebook we'll be exploring the urban heat island effect by looking at the impact on surface temperature of roof color and street trees. We'll be replicating the process described [here](https://urbanspatialanalysis.com/urban-heat-islands-street-trees-in-philadelphia/) but using Python tools rather than ESRI.\n",
     "\n",
-    "\n",
-    "**Extra packages:** To run this notebook, you'll need the PyViz tools and a library of *top of atmosphere* calculations from [`rio-toa`](https://github.com/mapbox/rio-toa): `pip install rio-toa`\n",
-    "\n",
-    "**Data sources:** This notebook uses Landsat data from Google Cloud Storage as well as some geographic data from OpenDataPhilly."
+    "**Data sources:** This notebook uses Landsat data from Google Cloud Storage as well as some geographic data from [OpenDataPhilly](https://opendataphilly.org/)."
    ]
   },
   {
@@ -26,33 +23,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import warnings\n",
-    "\n",
     "import intake\n",
     "import xarray as xr\n",
     "import pandas as pd\n",
     "import numpy as np\n",
-    "\n",
     "import geopandas as gpd\n",
-    "\n",
     "import cartopy.crs as ccrs\n",
-    "\n",
     "import hvplot.xarray  # noqa\n",
     "import hvplot.pandas  # noqa\n",
     "\n",
     "from geoviews.tile_sources import EsriImagery\n",
     "from pyproj import CRS"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "warnings.simplefilter('ignore')"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -246,7 +229,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds.crs"
+    "ds.rio.crs"
    ]
   },
   {
@@ -299,24 +282,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In this case, though, we are interested in a subset of data that covers that city of Philadelphia. So we need some geometry to specify the bounds of the city. We can get a GeoJSON of neighborhood data from [OpenDataPhilly](https://www.opendataphilly.org/dataset/philadelphia-neighborhoods/resource/6c61f240-aafe-478e-b993-b75fd09a93d6)."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "url = 'https://github.com/azavea/geo-data/raw/master/Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.geojson'\n",
-    "geoms = gpd.read_file(url)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We can compute the bounds of each of these neighborhoods and then using min and max get a rectangle that encompasses all of Philly."
+    "In this case, though, we are interested in a subset of data that covers that city of Philadelphia. So we need some geometry to specify the bounds of the city. We can get a GeoJSON of neighborhood data from [OpenDataPhilly](https://opendataphilly.org/datasets/philadelphia-neighborhoods/)."
    ]
   },
   {
@@ -325,6 +291,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "geoms = gpd.GeoDataFrame.from_features(cat.neighborhoods_philadelphia.read(), crs=4326)\n",
+    "\n",
     "bounds = geoms.geometry.bounds\n",
     "lower_left_corner_lat_lon = bounds.minx.min(), bounds.miny.min()\n",
     "upper_right_corner_lat_lon = bounds.maxx.max(), bounds.maxy.max()\n",
@@ -336,7 +304,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Using the crs defined above, we can transform these lat lons into map coordinates."
+    "Using the crs defined above, we can transform these lat/lons into map coordinates."
    ]
   },
   {
@@ -397,7 +365,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "band_plot = subset.mean('band').hvplot(x='x', y='y', datashade=True, project=True, crs=crs, cmap='gray')\n",
+    "geoms.head()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "band_plot = subset.mean('band').hvplot(x='x', y='y', rasterize=True, crs=crs, cnorm='eq_hist', cmap='gray')\n",
     "hood_plot = geoms.hvplot(geo=True, alpha=.5, c='mapname', legend=False, frame_height=450)\n",
     "\n",
     "band_plot * hood_plot"
@@ -428,7 +407,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In order to visualize NDVI, the data will need to be loaded and the NDVI computed. We can expect this to take some non-trivial amount of time (on the order of 20 sec on my machine)."
+    "In order to visualize NDVI, the data will need to be loaded and the NDVI computed."
    ]
   },
   {
@@ -437,7 +416,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "NDVI.hvplot(x='x', y='y', crs=crs, datashade=True, project=True, cmap='viridis', frame_height=450)"
+    "NDVI.hvplot(x='x', y='y', crs=crs, rasterize=True, cmap='viridis', frame_height=450)"
    ]
   },
   {
@@ -446,7 +425,7 @@
    "source": [
     "## Calculate land surface temperature\n",
     "\n",
-    "Given the NDVI calculated above, we can determine land surface temperature. For ease, we'll use some *top of atmosphere* calculations that have already been written up as Python functions as part of rasterio work in the `rio_toa` module. We'll also need to specify one more for transforming satellite temperature (brightness temp) to land surface temperature. For these calculations we'll use both Thermal Infrared bands - 10 and 11."
+    "Given the NDVI calculated above, we can determine land surface temperature. For ease, we'll use some *top of atmosphere* calculations that have already been written up as Python functions as part of rasterio work in the [rio_toa](https://github.com/mapbox/rio-toa/tree/master) module. We will reproduce the `brightness_temp` and `temp_rescale` functions from the `rio_toa` module here to make it easier and faster."
    ]
   },
   {
@@ -455,7 +434,55 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from rio_toa import brightness_temp, toa_utils"
+    "def brightness_temp(img, ML, AL, K1, K2, src_nodata=0):\n",
+    "    \"\"\"Calculate brightness temperature of Landsat 8 combining the calculations of radiance and temperature.\n",
+    "\n",
+    "    Parameters:\n",
+    "        img (ndarray): Input pixels array.\n",
+    "        ML (float): Multiplicative rescaling factor from scene metadata.\n",
+    "        AL (float): Additive rescaling factor from scene metadata.\n",
+    "        K1 (float): Thermal conversion constant from scene metadata.\n",
+    "        K2 (float): Thermal conversion constant from scene metadata.\n",
+    "        src_nodata (int, optional): Value representing 'no data' in input array.\n",
+    "\n",
+    "    Returns:\n",
+    "        ndarray: Array of at-satellite brightness temperature in degrees Kelvin.\n",
+    "    \"\"\"\n",
+    "    L = ML * img.astype(np.float32) + AL\n",
+    "    L[img == src_nodata] = np.nan\n",
+    "\n",
+    "    T = K2 / np.log((K1 / L) + 1)\n",
+    "\n",
+    "    return T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def temp_rescale(arr, temp_scale):\n",
+    "    \"\"\"Converts the `arr` in Kelvin (K) to the specified `temp_scale`- Kelvin(K), Farenheit(F) or Celcius(C).\"\"\"\n",
+    "    if temp_scale == 'K':\n",
+    "        return arr\n",
+    "\n",
+    "    elif temp_scale == 'F':\n",
+    "        return arr * (9 / 5.0) - 459.67\n",
+    "\n",
+    "    elif temp_scale == 'C':\n",
+    "        return arr - 273.15\n",
+    "\n",
+    "    else:\n",
+    "        raise ValueError('%s is not a valid temperature scale'\n",
+    "                         % (temp_scale))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We'll also need to specify one more for transforming satellite temperature (brightness temp) to land surface temperature. For these calculations we'll use both Thermal Infrared bands - 10 and 11."
    ]
   },
   {
@@ -513,10 +540,10 @@
     "    K1 = metadata[f'K1_CONSTANT_BAND_{band}']\n",
     "    K2 = metadata[f'K2_CONSTANT_BAND_{band}']\n",
     "    \n",
-    "    at_satellite_temp = brightness_temp.brightness_temp(data.sel(band=band).values, ML, AL, K1, K2)\n",
+    "    at_satellite_temp = brightness_temp(data.sel(band=band).values, ML, AL, K1, K2)\n",
     "    \n",
     "    temp = land_surface_temp(at_satellite_temp, w, NDVI).compute()\n",
-    "    return toa_utils.temp_rescale(temp, units)"
+    "    return temp_rescale(temp, units)"
    ]
   },
   {
@@ -592,12 +619,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "thresholded_temp_p = (mean_temp_f.where(mean_temp_f > 90)\n",
-    "    .hvplot(x='x', y='y', title='Mean Temp (F) > 90',\n",
-    "            cmap=special_cmap, crs=crs, frame_width=400,\n",
-    "            frame_height=450,\n",
-    "            colorbar=False, legend=False)\n",
-    "    .redim(value='Temperature (F)'))\n",
+    "thresholded_temp = mean_temp_f.where(mean_temp_f > 90)\n",
+    "    \n",
+    "thresholded_temp_p = thresholded_temp.hvplot(x='x', y='y', title='Mean Temp (F) > 90',\n",
+    "            cmap=special_cmap, crs=crs, frame_width=400, frame_height=450, legend=False)\n",
     "\n",
     "thresholded_temp_p + thresholded_temp_p.opts(alpha=.3, data_aspect=None) * EsriImagery"
    ]
@@ -616,43 +641,35 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
-    "# url = 'http://data.phl.opendata.arcgis.com/datasets/957f032f9c874327a1ad800abd887d17_0.geojson'\n",
-    "url = 'https://opendata.arcgis.com/api/v3/datasets/5abe042f2927486891c049cf064338cb_0/downloads/data?format=geojson&spatialRefId=4326&where=1%3D1'\n",
-    "trees_gdf = gpd.read_file(url)\n",
+    "trees_gdf = gpd.GeoDataFrame.from_features(cat.ppr_tree_inventory_2022.read(), crs=4326)\n",
     "trees_gdf.head()"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The trees plot can be generated straight from the geopandas dataframe, but that is rather slow. By inspecting the dataframe, we can see that each tree is represented as a point. We can create a simpler pandas dataframe with lat as one column and lon as the other."
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "trees_df = pd.DataFrame({'Longitude': trees_gdf.geometry.x, 'Latitude': trees_gdf.geometry.y})\n",
-    "trees_df.head()"
+    "tree_p = trees_gdf.hvplot.points('LOC_X', 'LOC_Y', title='Street Tree Density',\n",
+    "                                geo=True, rasterize=True, dynspread=True, max_px=5,\n",
+    "                                 threshold=0.9, frame_height=450, cmap='Greens', cnorm='eq_hist')\n",
+    "\n",
+    "thresholded_temp_p.opts(alpha=.5) * tree_p.opts(alpha=.5)"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "tree_p = trees_df.hvplot.points('Longitude', 'Latitude', title='Street Tree Density',\n",
-    "                                geo=True, datashade=True, dynspread=True, \n",
-    "                                frame_height=450)\n",
+    "Here, we use the `dynspread` parameter that helps us to dynamically calculate the points to display based on the local density of the points. This is more evident when you zoom into a sparse region of the plot. We also set the `threshold` parameter to control the fraction of the maximum density at which spreading starts, as well as `max_px` which sets the maximum number of pixels by which points can be spread.\n",
     "\n",
-    "thresholded_temp_p.opts(alpha=.5) * tree_p.opts(alpha=.5)"
+    "See [datashader.spreading](https://datashader.org/getting_started/Pipeline.html#spreading) and [dynspread](https://datashader.org/api.html#datashader.transfer_functions.dynspread) for more details."
    ]
   }
  ],
@@ -672,9 +689,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.18"
+   "version": "3.11.9"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }