Skip to content

Instantly share code, notes, and snippets.

@jdbcode

jdbcode/about_ee_with_drought_et.ipynb

Created April 14, 2024 14:06
Show Gist options
  • Save jdbcode/22558d59e9236da9e3e766e43a2007b6 to your computer and use it in GitHub Desktop.
Save jdbcode/22558d59e9236da9e3e766e43a2007b6 to your computer and use it in GitHub Desktop.
Download ZIP
about_ee_with_drought_et.ipynb
Display the source blob
Display the rendered blob
Raw
about_ee_with_drought_et.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"provenance": [],
"authorship_tag": "ABX9TyO0N/rGJf9h5D9eemjaebQC",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/jdbcode/22558d59e9236da9e3e766e43a2007b6/about_ee_with_drought_et.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"source": [
"#@title Welcome to Earth Engine! { display-mode: \"form\" }\n",
"%%HTML\n",
"\n",
"<div style=\"max-width: 960px;\">\n",
" <p style=\"font-size: 28px;\">\n",
" Google's mission is to <strong><span style=\"color: #4285F4;\">organize</span></strong> the world's <strong><span style=\"color: #EA4335;\">information</span></strong> and make it <strong><span style=\"color: #34A853;\">universally accessible</span></strong> and <strong><span style=\"color: #FBBC04;\">useful</span></strong>\n",
" </p>\n",
"\n",
" <img src=\"https://lh3.googleusercontent.com/u5uhPNPusadrV3H-OmU1KzOV3BTJaU82x9FDapFtkBFcFcNFj8Adu54VD0F6tvX7gh_M2_A59w7BIiY4m33ejA\" width=\"200\" style=\"float: left;\">\n",
"\n",
" <p style=\"font-size: 32px; margin: 40px 0px 0px 120px;\">\n",
" <strong><span style=\"color: #4285F4;\">Earth Engine</span></strong> is one of Google's platforms for organizing the world's <strong><span style=\"color: #34A853;\">environmental geospatial information</span></strong> and making it universally accessible and useful.\n",
" </p>\n",
"</div>"
],
"metadata": {
"id": "OkPyEvUKy2AB"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"#@title Earth Engine unlocks a massive information gateway { display-mode: \"form\" }\n",
"\n",
"import ipywidgets\n",
"from ipywidgets import VBox, HBox, Output, HTML\n",
"\n",
"video = HTML('''<iframe width=\"1000px\" height=\"564px\" src=\"https://www.youtube.com/embed/mlYxBLvLNzE?loop=1&autoplay=1&rel=0\" title=\"AMACRO, Brazil - Earth Timelapse\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture; web-share\" referrerpolicy=\"strict-origin-when-cross-origin\" allowfullscreen></iframe>''')\n",
"desc = HTML('''<h2 style=\"max-width: 500px; margin: 10px\">\n",
" With its petabyte-scale geospatial archive, powerful analysis tools, and\n",
" intuitive interface, Google Earth Engine enables us to monitor environmental\n",
" change, assess global trends, and make data-driven decisions that shape the\n",
" future of our planet.</h2><br>\n",
" <h2 style=\"max-width: 500px; margin: 10px\">Here we see 38 years of\n",
" deforestation in Brazil from Landsat imagery that was processed using\n",
" Earth Engine.</h2>''')\n",
"hbox = HBox([video, desc])\n",
"vbox = VBox([hbox])\n",
"vbox"
],
"metadata": {
"id": "28TX3Sdc-Ajf"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"# Earth Engine + Python\n",
"\n",
"### Earth Engine's Python client library brings the power of planetary-scale geospatial analysis directly to your fingertips. Harness Earth Engine's vast data and computation effortlessly within the familiar and versatile Python environment.\n",
"\n",
"### Earth Engine is installed by default in Google Colab – Google's Jupyter notebook environment – so all we have to do is import the Earth Engine library and connect to a Cloud project"
],
"metadata": {
"id": "twnN3iiJuUOb"
}
},
{
"cell_type": "code",
"source": [
"import ee\n",
"ee.Authenticate()\n",
"ee.Initialize(project='my-project')"
],
"metadata": {
"id": "UHh5Jd2Vuehc"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"#@title The **geemap** interface provides interactive data visualization and exploration { display-mode: \"form\" }\n",
"\n",
"import geemap\n",
"import geemap.colormaps as cm\n",
"\n",
"# Computation\n",
"m = geemap.Map(basemap='CartoDB.DarkMatter')\n",
"m.layout.max_width = '1000px'\n",
"\n",
"img = ee.ImageCollection(\"MODIS/061/MCD43A4\").filterDate('2018-04-01', '2018-06-01')\n",
"\n",
"trueColorVis = {\n",
" 'bands': ['Nadir_Reflectance_Band1', 'Nadir_Reflectance_Band4', 'Nadir_Reflectance_Band3'],\n",
" 'min': 0.0,\n",
" 'max': 4000.0,\n",
" 'gamma': 1.4,\n",
"}\n",
"\n",
"m.add_layer(img, trueColorVis, 'MODIS')\n",
"\n",
"\n",
"# Display\n",
"title = ipywidgets.HTML('<h1>Global imagery from NASA MODIS</h1>')\n",
"desc = ipywidgets.HTML('''<h2 style=\"max-width: 500px; margin: 10px\">\n",
" geemap uses the ipyleaflet library to display map tiles requested and served\n",
" from Earth Engine on demand. You can zoom and pan the map, change base\n",
" layers, edit map layer display properties, and inspect pixel values.</h2>''')\n",
"hbox = HBox([m, desc])\n",
"vbox = VBox([title, hbox])\n",
"vbox"
],
"metadata": {
"id": "cPDeGYSAu606"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"#@title Access and analyze over 1000 geospatial datasets from hundreds of providers { display-mode: \"form\" }\n",
"\n",
"stac = ipywidgets.HTML('''<iframe width=\"1000\" height=\"600\" src=\"https://radiantearth.github.io/stac-browser/#/external/storage.googleapis.com/earthengine-stac/catalog/catalog.json\" allow=\"https://colab.sandbox.google.com/*\"></iframe>''')\n",
"desc = ipywidgets.HTML('''<h2 style=\"max-width: 500px; margin: 10px\">\n",
" Earth Engine offers datasets spanning weather, climate, topography,\n",
" agriculture, landcover, and imagery, provided by government agencies like\n",
" NASA, international organizations like the World Resources Institute, and\n",
" research institutions such as the University of Maryland.</h2>\n",
" <h2 style=\"max-width: 500px; margin: 10px\">All of the data are organized\n",
" into a consistent format, optimized for Google's cloud infrastructure, and\n",
" colocated with compute servers to simplify your workflows and accelerates\n",
" analysis.</h2>''')\n",
"hbox = HBox([stac, desc])\n",
"vbox = VBox([hbox])\n",
"vbox"
],
"metadata": {
"id": "WjhdTQPQCwzm"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"#@title For instance, we can look at average drought conditions for 2023 { display-mode: \"form\" }\n",
"\n",
"# Computation\n",
"mean_pdsi = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE').filter(\n",
" ee.Filter.date('2023-01-01', '2024-01-01')).select('pdsi').mean()\n",
"\n",
"m = geemap.Map(basemap='CartoDB.DarkMatter', draw_ctrl=False, fullscreen_ctrl=False)\n",
"m.layout.max_width = '1000px'\n",
"vis_params = {\n",
" 'min': -1000,\n",
" 'max': 1000,\n",
" 'palette': cm.get_palette('RdBu', n_class=8),\n",
"}\n",
"m.add_layer(\n",
" mean_pdsi, vis_params, 'PDSI')\n",
"m.add_colorbar(\n",
" label='Drought index',\n",
" vis_params=vis_params)\n",
"\n",
"\n",
"# Display\n",
"title = ipywidgets.HTML('<h1>Mean Palmer Drought Severity Index (PDSI), 2023</h1>')\n",
"desc = ipywidgets.HTML('''<h2 style=\"max-width: 500px; margin: 10px\">\n",
" In 2023, stark global disparities in wet and dry conditions were evident.\n",
" Regions of Africa and South America experienced severe drought, underscoring\n",
" the critical need to monitor these environmental trends.</h2>''')\n",
"hbox = HBox([m, desc])\n",
"vbox = VBox([title, hbox])\n",
"vbox"
],
"metadata": {
"id": "UJS6GezOBC9g"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"#@title Python's matplotlib visualization tools make it easy to compare data across time { display-mode: \"form\" }\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.cm as cmap\n",
"import numpy as np\n",
"from mpl_toolkits.axes_grid1 import ImageGrid\n",
"\n",
"# Computation\n",
"years = ee.List.sequence(2020, 2023)\n",
"\n",
"def composite(year):\n",
" start = ee.Date.fromYMD(year, 1, 1)\n",
" end = ee.Date.fromYMD(year, 12, 31)\n",
" return ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE').filterDate(\n",
" start, end).select('pdsi').mean()\n",
"\n",
"band_names = years.getInfo()\n",
"band_names = [str(year) for year in band_names]\n",
"year_comps = ee.ImageCollection(years.map(composite)).toBands().rename(band_names)\n",
"\n",
"conus_mask = ee.ImageCollection(\"GRIDMET/DROUGHT\").first().select([0]).mask()\n",
"annual_drought = year_comps.updateMask(conus_mask).unmask(-9999).reproject('EPSG:5070', None, 10000).clipToBoundsAndScale(\n",
" geometry=ee.Geometry.BBox(-120, 25, -73, 50), scale=10000\n",
")\n",
"\n",
"annual_drought_npy = ee.data.computePixels({\n",
" 'expression': annual_drought,\n",
" 'fileFormat': 'NUMPY_NDARRAY',\n",
"})\n",
"\n",
"fig = plt.figure(figsize=(15.0, 15.0))\n",
"grid = ImageGrid(\n",
" fig,\n",
" 111,\n",
" nrows_ncols=(2, 2),\n",
" axes_pad=0.4,\n",
" cbar_mode=\"single\",\n",
" cbar_location=\"right\",\n",
" cbar_pad=0.4,\n",
" cbar_size=\"1.3%\",\n",
")\n",
"\n",
"for ax, name in zip(grid, band_names):\n",
" arr = annual_drought_npy[name]\n",
" arr = np.where(arr == -9999, np.nan, arr)\n",
" mask = np.isnan(arr)\n",
" masked_arr = np.ma.array(arr, mask=mask)\n",
" ax.imshow(masked_arr, vmin=-500, vmax=500, cmap=cmap.RdBu)\n",
" ax.set_facecolor('grey')\n",
" ax.set_title(name)\n",
"\n",
"colorbar = plt.colorbar(ax.get_children()[0], cax=grid[0].cax)\n",
"colorbar.set_label(\"Drought severity (unitless anomaly)\")\n",
"\n",
"\n",
"# Display\n",
"display(HTML('''<h1>PDSI 2020 to 2023 for the contiguous United States</h1>\n",
" <h2>With this plot we can compare drought conditions over time.</h2>'''))\n",
"plt.show()"
],
"metadata": {
"id": "17pbwsFc_4Zv"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"#@title Easily filter, query, and analyze data { display-mode: \"form\" }\n",
"\n",
"import pandas as pd\n",
"import altair as alt\n",
"\n",
"# Computation\n",
"img_col = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE').filter(\n",
" ee.Filter.date('2016-01-01', '2024-05-01')).select('pdsi')\n",
"aoi = ee.FeatureCollection('FAO/GAUL_SIMPLIFIED_500m/2015/level2').filter(\n",
" 'ADM2_CODE == 31640').geometry()\n",
"washington = ee.FeatureCollection('FAO/GAUL_SIMPLIFIED_500m/2015/level1').filter(\n",
" 'ADM1_CODE == 3261')\n",
"\n",
"def get_stats(geometry):\n",
" def calc(img):\n",
" stat = img.rename('PDSI').reduceRegion(\n",
" reducer=ee.Reducer.mean(),\n",
" geometry=geometry,\n",
" scale=10000)\n",
" return ee.Feature(None, stat).set(\n",
" {'Date': img.date().format('YYYY-MM-dd')})\n",
" return calc\n",
"\n",
"stats_fc = ee.FeatureCollection(img_col.map(get_stats(aoi)))\n",
"stats_fc\n",
"\n",
"pdsi_df = ee.data.computeFeatures({\n",
" 'expression': stats_fc,\n",
" 'fileFormat': 'PANDAS_DATAFRAME'\n",
"})\n",
"pdsi_df.drop('geo', axis=1, inplace=True)\n",
"\n",
"m = geemap.Map(lite_mode=True)\n",
"m.add_layer(aoi, {'color': '2166ac'})\n",
"m.layout.max_width = '420px'\n",
"m.layout.max_height = '400px'\n",
"m.center_object(washington, 6)\n",
"\n",
"\n",
"drought_chart = alt.Chart(pdsi_df).mark_bar().encode(\n",
" x=alt.X('Date:T'),\n",
" y=alt.Y('PDSI:Q'),\n",
" color=alt.Color(\n",
" 'PDSI:Q', scale=alt.Scale(scheme='redblue', domain=(-500, 500))),\n",
" tooltip=[\n",
" alt.Tooltip('Date:T'),\n",
" alt.Tooltip('PDSI:Q')\n",
" ]).properties(width=580, height=300)\n",
"\n",
"\n",
"# Display\n",
"table_output = Output()\n",
"chart_output = Output()\n",
"title = HTML('''<h1>What is the recent drought history for Franklin Country, WA?</h1>''')\n",
"desc = HTML('''<h2 style=\"max-width: 400px; margin: 20px 0px 0px 20px\">We can\n",
" filter drought data in Earth Engine by region and time and transfer the\n",
" results into a Pandas DataFrame and plot with Altair.<br><br>Notice the\n",
" significant drought in 2020 and 2021?</h2>''')\n",
"hbox = HBox([m, table_output, desc])\n",
"vbox = VBox([title, hbox, chart_output])\n",
"with table_output:\n",
" display(pdsi_df)\n",
"with chart_output:\n",
" display(drought_chart)\n",
"vbox"
],
"metadata": {
"id": "sz36P1mD4yIE"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"#@title Hone your analyses with Earth Engine's pixel-level algorithms { display-mode: \"form\" }\n",
"\n",
"# Computation\n",
"et = ee.ImageCollection(\n",
" 'OpenET/ENSEMBLE/CONUS/GRIDMET/MONTHLY/v2_0').select('et_ensemble_mad')\n",
"\n",
"years = ee.List.sequence(2016, 2023)\n",
"\n",
"def composite(year):\n",
" start = ee.Date.fromYMD(year, 1, 1)\n",
" end = ee.Date.fromYMD(ee.Number(year).add(1), 1, 1)\n",
" return ee.ImageCollection(\n",
" 'OpenET/ENSEMBLE/CONUS/GRIDMET/MONTHLY/v2_0').filterDate(\n",
" start, end).select('et_ensemble_mad').max().set('year', year)\n",
"\n",
"aoi = ee.FeatureCollection(\n",
" 'FAO/GAUL_SIMPLIFIED_500m/2015/level2').filter('ADM2_CODE == 31640').geometry()\n",
"\n",
"year_comps = ee.ImageCollection(years.map(composite))\n",
"\n",
"mean_et = year_comps.reduce(ee.Reducer.mean()).clip(aoi)\n",
"stdev_et = year_comps.reduce(ee.Reducer.stdDev()).clip(aoi)\n",
"\n",
"m = geemap.Map(basemap='OpenStreetMap.BlackAndWhite', center=[46.5, -118.8], zoom=9, lite_mode=True)\n",
"m.layout.max_width = '600px'\n",
"m.layout.max_height = '470px'\n",
"et_vis = {'min': 31.585, 'max': 215.665, 'palette': cm.get_palette('YlGnBu', n_class=8)}\n",
"m.add_layer(mean_et, et_vis, 'Mean of max annual ET (mm)')\n",
"m.add_colorbar(\n",
" label='Mean of max annual ET (mm)',\n",
" vis_params=et_vis)\n",
"\n",
"mean_et_hist = mean_et.reduceRegion(\n",
" reducer=ee.Reducer.autoHistogram(), geometry=aoi, scale=250)\n",
"data = mean_et_hist.getInfo()['et_ensemble_mad_mean']\n",
"df = pd.DataFrame(data, columns=['Mean of max annual ET (mm)', 'Regional frequency'])\n",
"\n",
"chart = alt.Chart(df).mark_bar().encode(\n",
" x=alt.X('Mean of max annual ET (mm):Q'),\n",
" y=alt.Y('Regional frequency:Q')\n",
")\n",
"\n",
"line_chart = alt.Chart().mark_rule(color='red').encode(\n",
" x=alt.X(datum=120)\n",
")\n",
"\n",
"not_irrigated = line_chart.mark_text(\n",
" x=\"width\",\n",
" dy=-120,\n",
" dx=-10,\n",
" align=\"right\",\n",
" baseline=\"bottom\",\n",
" text=\"Not irrigated\"\n",
")\n",
"\n",
"irrigated = line_chart.mark_text(\n",
" x=\"width\",\n",
" dy=-120,\n",
" dx=50,\n",
" align=\"right\",\n",
" baseline=\"bottom\",\n",
" text=\"Irrigated\"\n",
")\n",
"\n",
"chart = (chart + line_chart + not_irrigated + irrigated).properties(width=500, height=300)\n",
"\n",
"irrigated_mask = mean_et.gt(120)\n",
"m1 = geemap.Map(basemap='OpenStreetMap.BlackAndWhite', center=[46.5, -118.8], zoom=9, lite_mode=True)\n",
"m1.layout.max_width = '1180px'\n",
"m1.layout.max_height = '370px'\n",
"m1.add_layer(aoi, {'color': 'white'}, 'Franklin County, WA')\n",
"m1.add_layer(mean_et.updateMask(irrigated_mask), et_vis, 'Mean of max annual ET (mm)')\n",
"\n",
"# Display\n",
"chart_output = Output()\n",
"title = HTML('''<h1>Using evapotranspiration (ET), can we separate irrigated land from non-irrigated land</h1>''')\n",
"hbox = HBox([m, VBox([chart_output, HTML('''<h2 style=\"max-width: 500px; margin: 10px\">\n",
" We can see from the map and histogram of mean max evapotranpiration in\n",
" Franklin County that irrigated and non-irrigated lands have very distinct\n",
" ranges of max ET. We can use the value of max ET 120 mm as a divider.</h2>''')])])\n",
"title2 = HTML('''<h2>Only irrigated lands in Franklin Country, WA</h2>''')\n",
"vbox = VBox(children=[title, hbox, title2, m1])\n",
"with chart_output:\n",
" display(chart)\n",
"display(vbox)"
],
"metadata": {
"id": "lolZckVjXIWb"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"#@title Gain insights needed to make informed, strategic decisions { display-mode: \"form\" }\n",
"\n",
"# Computation\n",
"drought_filter = ee.Filter.inList(leftField='year', rightValue=[2020, 2021])\n",
"drought_years = year_comps.filter(drought_filter).mean()\n",
"not_drought_years = year_comps.filter(drought_filter.Not()).mean()\n",
"\n",
"dif = drought_years.subtract(not_drought_years)\n",
"dif_rel = drought_years.divide(not_drought_years).multiply(100).updateMask(irrigated_mask).clip(aoi)\n",
"\n",
"m = geemap.Map(basemap='OpenStreetMap.BlackAndWhite', center=[46.5, -118.8], zoom=9, lite_mode=True)\n",
"m.layout.max_width = '1180px'\n",
"m.layout.max_height = '370px'\n",
"m.add_layer(aoi, {'color': 'white'}, 'Franklin County, WA')\n",
"et_dif_vis = {'min': 80, 'max': 120, 'palette': cm.get_palette('RdBu', n_class=8)}\n",
"m.add_layer(dif_rel, et_dif_vis, 'Drought ET relative to non-drought (%)')\n",
"m.add_colorbar(\n",
" label='Drought ET relative to non-drought (%)',\n",
" vis_params=et_dif_vis)\n",
"\n",
"mean_et_hist = dif_rel.reduceRegion(reducer=ee.Reducer.autoHistogram(), geometry=aoi, scale=100)\n",
"data = mean_et_hist.getInfo()['et_ensemble_mad']\n",
"df = pd.DataFrame(data, columns=['x', 'y'])\n",
"\n",
"chart = alt.Chart(df).mark_bar(size=8).encode(\n",
" x=alt.X('x:Q',\n",
" title='Drought period mean max ET relative to non-drought (%)',\n",
" scale=alt.Scale(domain=[75, 130], clamp=True)),\n",
" y=alt.Y('y:Q',\n",
" title='Regional frequency'),\n",
" color=alt.Color('x',\n",
" title='%',\n",
" scale=alt.Scale(scheme='redblue', domain=[80, 120]))\n",
").properties(width=600, height=300)\n",
"\n",
"line_chart = alt.Chart().mark_rule(color='black').encode(\n",
" x=alt.X(datum=100)\n",
")\n",
"\n",
"less_et = line_chart.mark_text(\n",
" x=\"width\",\n",
" dy=-130,\n",
" dx=-10,\n",
" align=\"right\",\n",
" baseline=\"bottom\",\n",
" text=\"Less ET\"\n",
")\n",
"\n",
"more_et = line_chart.mark_text(\n",
" x=\"width\",\n",
" dy=-130,\n",
" dx=50,\n",
" align=\"right\",\n",
" baseline=\"bottom\",\n",
" text=\"More ET\"\n",
")\n",
"\n",
"final_chart = chart + line_chart + less_et + more_et\n",
"\n",
"# Display\n",
"output1 = Output()\n",
"output2 = Output()\n",
"desc = HTML('''<h1 style=\"max-width: 1000px\">How does drought\n",
" affect max evapotranspiration for irrigated lands?<h1>\n",
" <h2 style=\"max-width: 1000px\">Evapotranspiration is\n",
" the loss of water from land and plants to the atmosphere. Understanding\n",
" how drought affects ET is important for balancing crop productivity and\n",
" water conservation. We can isolate ET for periods of drought and not and\n",
" compare them. Here, we calculate the mean max ET for drought and non-drought\n",
" periods and compare the results by dividing drought ET by non-drought ET.\n",
" The result indicates how different drought max ET is from non-drought\n",
" max ET.</h2>''')\n",
"hbox = HBox([output1, output2])\n",
"vbox = VBox(children=[desc, m, hbox])\n",
"display(vbox)\n",
"with output1:\n",
" display(final_chart)\n",
"\n",
"with output2:\n",
" display(HTML('''<h2 style=\"max-width: 380px; margin: 10px\">The majority of\n",
" irrigated lands in Franklin County, WA have greater or equal max\n",
" evapotranspiration during drought periods relative to non-drought periods</h2>'''))"
],
"metadata": {
"id": "VNiflmf00-kR"
},
"execution_count": null,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment