diff --git a/notebooks/heat_islands/zBackup/Exploring_Urban_Heat_with_ECOSTRESS_LST.ipynb b/notebooks/heat_islands/Exploring_Urban_Heat_with_ECOSTRESS_LST.ipynb similarity index 56% rename from notebooks/heat_islands/zBackup/Exploring_Urban_Heat_with_ECOSTRESS_LST.ipynb rename to notebooks/heat_islands/Exploring_Urban_Heat_with_ECOSTRESS_LST.ipynb index 789e4fb..64e8401 100644 --- a/notebooks/heat_islands/zBackup/Exploring_Urban_Heat_with_ECOSTRESS_LST.ipynb +++ b/notebooks/heat_islands/Exploring_Urban_Heat_with_ECOSTRESS_LST.ipynb @@ -7,7 +7,7 @@ "source": [ "## Application of ECOSTRESS Data: Exploring the imapcts of Urban Heat Islands\n", "\n", - "ADD NOTES HERE" + "In this notebook, we will demonstrate a simple example of using ECOSTRESS data together with GIS layers for the Boulder Colorado area to investigate spatial variation in Land Surface Temperature (LST). We will explore temperature variation within the Boulder, CO area as well as explore the relationship of LST and tree cover using the Boulder CO Open Tree Data" ] }, { @@ -31,13 +31,18 @@ "import os, time, shutil\n", "import warnings\n", "\n", + "import numpy as np\n", "import cartopy.crs as ccrs\n", "import rioxarray as rxr\n", "import hvplot.xarray\n", "import hvplot.pandas\n", "import holoviews as hv\n", + "import matplotlib.pyplot as plt\n", + "import rasterio\n", + "import rasterio.plot\n", "import geoviews as gv\n", "import geopandas as gpd\n", + "import shapely\n", "\n", "# Some cells may generate warnings that we can ignore. Comment below lines to see.\n", "warnings.filterwarnings('ignore')\n", @@ -91,7 +96,7 @@ "source": [ "### Step 2. Load Data\n", "\n", - "ADD NOTES HERE" + "We will now load the notebook data, including the GIS layers and the ECOSTRESS image" ] }, { @@ -127,17 +132,22 @@ "source": [ "### Load the tree inventory data\n", "# Define a path to the open tree data\n", - "tree_path = os.path.join(os.path.expanduser(\"~\"),\"HYR-SENSE\",\n", + "trees_path = os.path.join(os.path.expanduser(\"~\"),\"HYR-SENSE\",\n", " # GIS directory\n", " 'data', 'Urban_heat',\n", " # City of Boulder file\n", " 'Tree_Inventory_Open_Data.zip'\n", ")\n", "# Read file and merge geometries\n", - "tree_gdf = gpd.read_file(tree_path).dissolve()\n", + "#trees_gdf = gpd.read_file(trees_path).dissolve()\n", + "trees_gdf = gpd.read_file(trees_path)\n", + "trees_gdf = trees_gdf[trees_gdf.geometry!=None]\n", + "\n", + "# Clip to the ROI\n", + "trees_gdf = trees_gdf.sjoin(boulder_gdf)\n", "\n", "# Check that we have a GeoDataFrame\n", - "tree_gdf" + "trees_gdf" ] }, { @@ -169,7 +179,7 @@ "id": "9b074cd3-e830-454e-b3c8-ae0df12efcf2", "metadata": {}, "source": [ - "### Step 3. Prepare the ECOSTRESS data for analysis\n", + "### Step 3. Explore the GIS layers\n", "\n", "ADD NOTES HERE" ] @@ -190,18 +200,13 @@ "outputs": [], "source": [ "### Plot the GIS boundary for Boulder, Colorado\n", - "(\n", - " boulder_gdf.hvplot(\n", - " **map_opts,\n", - " line_color='orange', line_width=3,\n", - " fill_color='white', fill_alpha=.05,\n", - " )\n", - " *\n", - " gpd.GeoDataFrame(geometry=boulder_gdf.envelope).hvplot(\n", - " **map_opts,\n", - " line_color='skyblue', fill_color=None\n", - " )\n", - ").opts(**size_opts)" + "boulder_roi_plot = (boulder_gdf.hvplot(**map_opts,\n", + " line_color='orange', line_width=3,\n", + " fill_color='white', fill_alpha=.05) * \n", + " gpd.GeoDataFrame(geometry=boulder_gdf.envelope).hvplot(**map_opts,\n", + " line_color='skyblue', fill_color=None)\n", + ").opts(**size_opts)\n", + "boulder_roi_plot" ] }, { @@ -212,13 +217,21 @@ "outputs": [], "source": [ "### Plot the tree locations\n", - "trees_plot = tree_gdf.hvplot.points(\n", + "trees_plot = trees_gdf.hvplot.points(\n", " **map_opts,\n", " size=15\n", ").opts(**size_opts)\n", "trees_plot" ] }, + { + "cell_type": "markdown", + "id": "81fc66ea-62d9-4c5e-9282-19ef0688b26c", + "metadata": {}, + "source": [ + "### Step 4. Prepare the ECOSTRESS data for analysis" + ] + }, { "cell_type": "markdown", "id": "789fe8ff-e9aa-48b2-a9ab-2648b9be4984", @@ -235,16 +248,13 @@ "outputs": [], "source": [ "### Load the ECOSTRESS data\n", - "(\n", - " eco_lst_ds\n", - " .rio.reproject('EPSG:4326')\n", - " .hvplot.image(\n", - " x='x', y='y', **size_opts, \n", - " cmap='inferno', tiles='ESRI', \n", - " xlabel='Longitude', ylabel='Latitude', \n", - " title='ECOSTRESS LST (K)', \n", - " crs='EPSG:4326')\n", - ")\n" + "eco_lst_ds_plot = (eco_lst_ds.rio.reproject('EPSG:4326')\n", + " .hvplot.image(x='x', y='y', **size_opts, \n", + " cmap='inferno', tiles='ESRI', \n", + " xlabel='Longitude', ylabel='Latitude', \n", + " title='ECOSTRESS LST (K)', \n", + " crs='EPSG:4326'))\n", + "eco_lst_ds_plot" ] }, { @@ -276,16 +286,24 @@ "outputs": [], "source": [ "### Confirm we have data in degrees C\n", - "(\n", - " eco_lst_ds_degC\n", - " .rio.reproject('EPSG:4326')\n", - " .hvplot.image(\n", - " x='x', y='y', **size_opts, \n", - " cmap='inferno', tiles='ESRI', \n", - " xlabel='Longitude', ylabel='Latitude', \n", - " title='ECOSTRESS LST (degC)', \n", - " crs='EPSG:4326')\n", - ")" + "eco_lst_ds_degC_plot = (eco_lst_ds_degC.rio.reproject('EPSG:4326')\n", + " .hvplot.image(x='x', y='y', **size_opts, \n", + " cmap='inferno', tiles='ESRI', \n", + " xlabel='Longitude', ylabel='Latitude', \n", + " title='ECOSTRESS LST (degC)', \n", + " crs='EPSG:4326'))\n", + "eco_lst_ds_degC_plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "79417bff-24de-49a6-a7a0-49148fe417ff", + "metadata": {}, + "outputs": [], + "source": [ + "# Cleanup - remove original ECOSTRESS LST object in degrees Kelvin\n", + "del eco_lst_ds" ] }, { @@ -368,7 +386,7 @@ " line_color='skyblue', fill_color=None\n", " )\n", ").opts(**size_opts)\n", - "trees_plot = tree_gdf.hvplot.points(\n", + "trees_plot = trees_gdf.hvplot.points(\n", " **map_opts,\n", " size=15\n", ").opts(**size_opts)\n", @@ -376,31 +394,162 @@ "combined_plot" ] }, + { + "cell_type": "markdown", + "id": "9da24cef-14e1-4003-8623-0f76e3f84607", + "metadata": {}, + "source": [ + "In the interactive plot above you can pan around the ECOSTRESS LST image as well as zoom into the Boulder Colorado area to view the LST variation across the city and see the location of trees provided by the Open Tree Inventory Data. What patterns can you see in the larger image moving from west to east across the image which reflects a change in elevation and increasing population density? What patterns can you see as you zoom into Boulder and in relation to the location of dense tree cover?" + ] + }, { "cell_type": "markdown", "id": "ed8f4bf8-83a7-43d9-81a7-a32895d85b5f", "metadata": {}, "source": [ - "### Step 4. Analyze the ECOSTRESS data" + "### Step 5. Analyze the ECOSTRESS data\n", + "\n", + "Next, we will dive deeper into the data. First let's generate some statistics using the GIS layers and the ECOSTRESS image" + ] + }, + { + "cell_type": "markdown", + "id": "8556a510-302d-40cd-82ae-5d6229fc2d65", + "metadata": {}, + "source": [ + "#### Step 5a. Explore LST variation across the image using transects" ] }, { "cell_type": "code", "execution_count": null, - "id": "ab0d5b26-6460-46e9-ae53-a08f7fde9110", + "id": "a8c9c375-7976-4bf0-a268-f5894159dc88", "metadata": {}, "outputs": [], "source": [ - "#https://py.geocompx.org/05-raster-vector" + "# Create a transect from west to east that intersects with the Boulder city ROI - we will use this to create an LST transect plot\n", + "coords = [[-106.0, 40.02], [-105.0, 40.02]]\n", + "lst_transect = shapely.LineString(coords)\n", + "lst_transect_utm = gpd.GeoSeries(lst_transect, crs=4326).to_crs(32613)\n", + "lst_transect_utm = lst_transect_utm.iloc[0]\n", + "distances = np.arange(0, lst_transect_utm.length, 250)\n", + "lst_transect_pnt = [lst_transect_utm.interpolate(distance) for distance in distances]\n", + "lst_transect_pnt = gpd.GeoSeries(lst_transect_pnt, crs=32613).to_crs('EPSG:4326')\n", + "# Convert to a geo data frame\n", + "lst_transect_pnt = gpd.GeoDataFrame(geometry=lst_transect_pnt)\n", + "lst_transect_pnt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e633937a-55a9-45f1-95d7-27ddb2ec96a0", + "metadata": {}, + "outputs": [], + "source": [ + "# Extract values of LST to each point along our previously generated transect of points\n", + "lst_transect_pnt['LST'] = lst_transect_pnt.geometry.apply(\n", + " lambda geom: eco_lst_ds_degC.rio.reproject('EPSG:4326').sel(x=geom.x, y=geom.y, \n", + " method=\"nearest\").values)\n", + "# Insert the distance data from west to east in meters by 250m steps\n", + "lst_transect_pnt['dist'] = distances\n", + "lst_transect_pnt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "521b9ce2-7789-4fdb-83a9-e50a68cbada3", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a figure showing the transect and the profile of LST values along the transect\n", + "fig, ax = plt.subplots(figsize=(8,6))\n", + "np.clip(eco_lst_ds_degC.rio.reproject('EPSG:4326').\n", + " squeeze(),0,30).plot.imshow(cmap='jet')\n", + "plt.title(\"ECOSTRESS LST Transect\")\n", + "gpd.GeoSeries(lst_transect).plot(ax=ax, color='black', \n", + " linewidth=2)\n", + "boulder_gdf.plot(ax=ax, color='none', edgecolor='white');\n", + "\n", + "# LST profile\n", + "fig, ax = plt.subplots(figsize=(8,6))\n", + "lst_transect_pnt.set_index('dist')['LST'].plot(ax=ax, \n", + " color='black',\n", + " linewidth=2)\n", + "ax.set_xlabel('Distance (m)')\n", + "ax.set_ylabel('LST (degC)');" + ] + }, + { + "cell_type": "markdown", + "id": "3c39f21a-9adc-4546-9579-0ec9552a8a3d", + "metadata": {}, + "source": [ + "#### Step 5b. Crop the ECOSTRESS data to the Boulder CO " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "081236b5-900f-45f5-862a-bb3bc842c7c5", + "metadata": {}, + "outputs": [], + "source": [ + "# Crop ECOSTRESS data to Region of Interest\n", + "# ---\n", + "\n", + "boulder_lst_da = (eco_lst_ds_degC\n", + " .rio.reproject('EPSG:4326')\n", + " .rio.clip(boulder_gdf.\n", + " to_crs('EPSG:4326').geometry))" ] }, { "cell_type": "code", "execution_count": null, - "id": "88ffd9a3-0c23-4070-9951-62f717f9fc55", + "id": "18846f0e-801a-493e-93c5-69af7fc56974", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "### Plot the data cropped to the Boulder CO ROI\n", + "eco_plot = (\n", + " boulder_lst_da\n", + " .rio.reproject('EPSG:4326')\n", + " .hvplot.image(\n", + " x='x', y='y', **size_opts, \n", + " cmap='inferno', tiles='ESRI', \n", + " xlabel='Longitude', ylabel='Latitude', \n", + " title='ECOSTRESS LST (degC)', \n", + " crs='EPSG:4326')\n", + ")\n", + "gis_boundary_plot = (\n", + " boulder_gdf.hvplot(\n", + " **map_opts,\n", + " line_color='black', line_width=5,\n", + " fill_color='white', fill_alpha=.05,\n", + " )\n", + " *\n", + " gpd.GeoDataFrame(geometry=boulder_gdf.envelope).hvplot(\n", + " **map_opts,\n", + " line_color='skyblue', fill_color=None\n", + " )\n", + ").opts(**size_opts)\n", + "trees_plot = trees_gdf.hvplot.points(\n", + " **map_opts,\n", + " size=3\n", + ").opts(**size_opts)\n", + "combined_plot = eco_plot * gis_boundary_plot * trees_plot\n", + "combined_plot" + ] + }, + { + "cell_type": "markdown", + "id": "a06e2735-9c22-4279-bd30-9bebe4a43cb9", + "metadata": {}, + "source": [ + "#### Step 5c. Summarize the LST data for the tree cover areas" + ] }, { "cell_type": "code", @@ -408,7 +557,59 @@ "id": "2c593675-6e1c-4653-9bd5-050a1a46f3d6", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "# Find the LST nearest each tree location\n", + "# ---\n", + "\n", + "trees_gdf['LST'] = trees_gdf.geometry.apply(\n", + " lambda geom: boulder_lst_da.sel(x=geom.x, y=geom.y, method=\"nearest\").values)\n", + "trees_gdf" + ] + }, + { + "cell_type": "markdown", + "id": "c851a88f-49f6-488f-88f9-886f23595e5b", + "metadata": {}, + "source": [ + "Review the table above which contains the information on individual trees in Boulder CO, their genus/species, and other information, and a new column added to the table containing the ECOSTRESS LST values in degrees C for the closest pixel to the location of each tree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb0425c0-ef8e-4e01-9fe4-2837273dd510", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot LST distribution for trees vs. total LST distribution\n", + "# ---\n", + "\n", + "tree_lst_plot = (trees_gdf.LST.hvplot.violin(frame_width=600, frame_height=400, fontscale=1.5,\n", + " label='', violin_fill_color='forestgreen', box_fill_color='forestgreen')\n", + " * boulder_lst_da.to_dataframe(name='LST').LST.hvplot.violin(\n", + " label='Total', violin_fill_color='skyblue', box_fill_color='skyblue'))\n", + "tree_lst_plot" + ] + }, + { + "cell_type": "markdown", + "id": "d4254e8d-3ff8-49fb-b1da-27121207c198", + "metadata": {}, + "source": [ + "### What do you observe?\n", + "\n", + "*Write your response here*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab0d5b26-6460-46e9-ae53-a08f7fde9110", + "metadata": {}, + "outputs": [], + "source": [ + "#https://py.geocompx.org/05-raster-vector" + ] } ], "metadata": {