Week1_ERA5Maps_GEE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dbd11bc3-d2f8-4257-8887-3268f0ecf2f4",
   "metadata": {},
   "source": [
    "# EDS 220 Fall 2021\n",
    "\n",
    "## Visualizing global precipitation using Google Earth Engine\n",
    "\n",
    "Let's get started with some data visualization and processing! Our tool for this exercise will be the Google Earth Engine (GEE):\n",
    "https://earthengine.google.com/\n",
    "\n",
    "GEE is a really powerful platform that allows easy access to many datasets without the need for a lot of pre-processing, so it's a good way to get our hands dirty with some basic mapping skills. Here we'll use the Python API, which will let us code in 'familiar' Python commands while still interfacing with the GEE platform. \n",
    "\n",
    "#### Recap from last time (GEE_Setup.ipynb)\n",
    "Prior to starting this notebook, you should have gone through the following steps:\n",
    "1. **Activating the Python environment containing your local GEE installation**\n",
    "   (on the command line in your terminal window) conda activate gee_env  \n",
    "   where \"gee_env\" is the name you gave to the environment creates for GEE to live in\n",
    "     \n",
    "2. **Authenticating Earth Engine**\n",
    "   (also on the terminal command line) earthengine authenticate  \n",
    "   This should bring up a login window in your Internet browser, where after entering your Google account information it will give you an authorization code to enter in the terminal\n",
    "     \n",
    "3. **Launching JupyterLab from within the appropriate Python environment**  \n",
    "   - If starting JupyterLab from the command line: make sure it says \"(gee_env)\" before your login name. This happens after the \"conda activate gee_env\" step, and your terminal prompt will look something like:  \n",
    "   (gee_env) Samanthas-MacBook-Pro:~ samstevenson$    \n",
    "     \n",
    "   - If starting JupyterLab from the Anaconda Navigator: make sure to go to the drop-down menu that says \"Applications on...\" and select \"gee_env\" before launching JupyterLab\n",
    "\n",
    "\n",
    "#### Importing Packages\n",
    "The code chunk below imports the needed packages: \n",
    "- **ee**, which covers the Google Earth Engine interface\n",
    "- **geemap**, which covers the GEE mapping functionality\n",
    "- **pandas**, to allow for manipulation of data subsets later on\n",
    "\n",
    "They're all included at the start of the notebook here, since it's good practice to do all of the importing of packages in one place (best at the beginning of the code) so you can easily see what tools are available to Python.\n",
    "\n",
    "As we discussed in the GEE_Setup notebook, here 'pandas' is being imported with a shorter name 'pd' to make it easier to reference that package later on. This is done with the \"import...as\" syntax. note that you *could* do this with any of the other packages as well, and you can call them whatever you want!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "36f523a9-c8a8-4130-bb5e-a186c73d7887",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import geemap\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b8e4852-1bd2-4c39-acd4-a5a8050e6389",
   "metadata": {},
   "source": [
    "## Create an interactive map\n",
    "\n",
    "Now let's load the package that will allow us to map a given dataset. The default basemap is, unsurprisingly, Google Maps. The code chunk below will display an empty Google Map that you can manipulate just like a regular Google Maps interface.  \n",
    "\n",
    "We accomplish this by calling the \"Map\" function within the geemap package: this uses the standard Python syntax of  \n",
    "[package name].[function name]  \n",
    "\n",
    "... or in this case,  \n",
    "`geemap.Map`  \n",
    "\n",
    "The output will be stored in an \"object\", which we will be able to do various things to in subsequent code chunks. (This is standard practice for \"object oriented languages\" like Python!)  \n",
    "\n",
    "Also important: below, we are providing an *argument* to the Map function. This is telling it to do two different things: \n",
    "- center the map at a particular lat/lon point: in this case, 40N and 100E\n",
    "- set the *zoom level* to a particular value. Large numbers mean more magnification (i.e. smaller fields of view). Here I set it to 2 since that seemed to be visually pretty appealing, but you can pick whatever you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fe2f5bfd-5a62-4ca9-b64e-0e5b4c16be81",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c8da3429df5f475383fe784c3bb72a67",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[40, -100], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(T…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Map = geemap.Map(center=[40,-100], zoom=2)\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8053919-9ad6-482b-a984-a48dc44e6a61",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4be51235-5b58-45d1-9a4e-2ba3b33b758f",
   "metadata": {},
   "source": [
    "## Load data and plot over the basemap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ffa288e-e325-446d-a147-97be7e5c5292",
   "metadata": {},
   "source": [
    "Now we need to tell GEE which dataset to layer on top of the base map. \n",
    "\n",
    "The next chunk of code will load some data; when executed, it will display within the map window above.\n",
    "\n",
    "#### ERA5 Reanalysis Data\n",
    "The dataset we'll be using here is the ERA5 \"reanalysis\": essentially, a reanalysis is a numerical model (in this case, a weather model) that's been forced to match observations of the real world (in this case, atmospheric weather observations) as closely as possible. ERA5 is produced by the European Centre for Medium-Range Weather Forecasts, or ECMWF; so you'll see below that the parent directory for this dataset is called ECMWF. A description of the dataset is here:\n",
    "https://developers.google.com/earth-engine/datasets/catalog/ECMWF_ERA5_DAILY\n",
    "\n",
    "##### Load in ERA5: GEE Image Collections\n",
    "The below chunk of code loads in the \"mean_2m_air_temperature\" and \"total precipitation\" fields, which are basically what they sound like: the temperature of air at an altitude of 2 meters, and the daily total precipitation, everywhere on Earth. Here we'll be working just with precipitation, but I've included temperature as well so you can see that other variables are just as easy to grab!\n",
    "\n",
    "There are a couple of things to know about how GEE stores image data. The ImageCollection function grabs a *set* of individual images that satisfies some criterion that you pass to GEE through the \"ee\" package; this is stored as an \"ImageCollection\" object which can be filtered and processed in various ways later.  \n",
    "\n",
    "The syntax for this is:  \n",
    "`ee.ImageCollection([args])`\n",
    "\n",
    "where you can pass different arguments to tell GEE which ImageCollection you would like. Here I'm going to tell it that I would like all of the ERA5 data: that means passing a location called \"ECMWF/ERA5/Daily\" that I got by poking through the [Earth Engine data catalog](https://developers.google.com/earth-engine/datasets/). There are also other ERA5 datasets: since we're going to explore precipitation, I'm using daily information so we can see individual rain events.\n",
    "\n",
    "*Note*: I'm also including comments within the code chunks themselves, noting what they do: this doesn't matter as much within a Jupyter notebook, but if you're writing a Python script, it's important to note what the various lines of code mean for your later self and anyone else who reads your code!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fe3c487e-4238-4a9e-a580-e02774c7e809",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add Earth Engine dataset: daily data from ERA5\n",
    "gdat = ee.ImageCollection('ECMWF/ERA5/DAILY')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c6b686e-c52a-4b7b-9ec1-4ce6e2c42e66",
   "metadata": {},
   "source": [
    "#### Select an image to plot\n",
    "\n",
    "To plot a map over the Google Maps basemap, we need not an \"ImageCollection\", but an \"Image\". The following code chunks go through the process of:\n",
    "- selecting a variable to plot (precipitation)\n",
    "- filtering that variable to a particular time range\n",
    "- calculating the average over that time range\n",
    "- layering the data on top of the basemap\n",
    "\n",
    "The first thing we need to do is to pick what quantity we would like to plot: ERA5 contains many different climate variables. Below are two examples of selecting a field from an ImageCollection: I'm including two so you can get a sense of how it might work, but we'll be focusing on precipitation since it's more visually interesting.\n",
    "\n",
    "The selection process is accomplished using the \"select\" function, applied to the \"gdat\" object we created in the code chunk above.\n",
    "\n",
    "**Recall:** the object \"gdat\" is a *collection* of images, that we will use the `.select` function to pick from:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "63d40aec-bcb6-4c68-be1e-ca2b8f63ce40",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select temperature and precipitation data from the ERA5 daily images\n",
    "ts = gdat.select('mean_2m_air_temperature')\n",
    "\n",
    "pr=gdat.select('total_precipitation')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49ce01f5-9b8d-4981-b654-31c98dac2459",
   "metadata": {},
   "source": [
    "Let's take a look at the \"pr\" object for a second. This can be done using the `print` command, or just by simply typing the name of the object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fc443a4d-c5d3-4a65-a9c8-b126f10ae09d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ee.ImageCollection({\n",
      "  \"functionInvocationValue\": {\n",
      "    \"functionName\": \"Collection.map\",\n",
      "    \"arguments\": {\n",
      "      \"baseAlgorithm\": {\n",
      "        \"functionDefinitionValue\": {\n",
      "          \"argumentNames\": [\n",
      "            \"_MAPPING_VAR_0_0\"\n",
      "          ],\n",
      "          \"body\": {\n",
      "            \"functionInvocationValue\": {\n",
      "              \"functionName\": \"Image.select\",\n",
      "              \"arguments\": {\n",
      "                \"bandSelectors\": {\n",
      "                  \"constantValue\": [\n",
      "                    \"total_precipitation\"\n",
      "                  ]\n",
      "                },\n",
      "                \"input\": {\n",
      "                  \"argumentReference\": \"_MAPPING_VAR_0_0\"\n",
      "                }\n",
      "              }\n",
      "            }\n",
      "          }\n",
      "        }\n",
      "      },\n",
      "      \"collection\": {\n",
      "        \"functionInvocationValue\": {\n",
      "          \"functionName\": \"ImageCollection.load\",\n",
      "          \"arguments\": {\n",
      "            \"id\": {\n",
      "              \"constantValue\": \"ECMWF/ERA5/DAILY\"\n",
      "            }\n",
      "          }\n",
      "        }\n",
      "      }\n",
      "    }\n",
      "  }\n",
      "})\n"
     ]
    }
   ],
   "source": [
    "pr\n",
    "print(pr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dee26d3-a4d1-42c5-b90a-a06a9d5cfe59",
   "metadata": {},
   "source": [
    "As you can see, `pr` is still an ImageCollection. This makes it not very mappable: in fact, this object contains information every day from 1979 to 2020. In order to create a map, we need to filter it down to a single field for a time of our choice. \n",
    "\n",
    "Let's select an arbitrary time range: here, December 1-2, 2019.\n",
    "\n",
    "The below code might look confusing to start with, so I'll break it down a bit: the first thing it does is applies the Python built-in `filter` command to our precipitation ImageCollection object. `filter` is a standard Python function: essentially, what it does is\n",
    "- take any object you give it (in this case, \"pr\")\n",
    "- apply any arbitrary function to that data\n",
    "\n",
    "Here, the function we're passing it is confusingly ALSO called filter, but it's the filtering function that's part of the \"ee\" package. This is `ee.Filter.date`, which takes a range of dates as arguments.\n",
    "\n",
    "**The point is** that the piece of the code that looks like this:\n",
    "`pr.filter(ee.Filter.date('2019-12-01', '2019-12-02'))`\n",
    "\n",
    "is saying \"take the pr object, and use Google Earth Engine's date filtering functionality to pick out the bits between December 1 and 2, 2019\".\n",
    "\n",
    "#### Other nifty Python thing: multiple functions applied in one line of code\n",
    "\n",
    "You'll notice that the line of code below ALSO includes another function: `.mean()`. This is exactly what it sounds like: an averaging operator. `.mean()` is another standard Python built-in, that takes whatever is before the period (in this case, precipitation filtered between December 1-2, 2019) and calculates the average.\n",
    "\n",
    "**The end result**: In one line of code, we can grab all of the precipitation data in a particular time range, and calculate the average, to get a single image that we can then map!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ecb38bd3-9745-4de6-a7c0-81ff94a56bf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter temperature and precipitation according to a specified time range, calculate the mean\n",
    "prdflt=pr.filter(ee.Filter.date('2019-12-01', '2019-12-02')).mean();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9343535-5832-4930-8122-a0ce38b0f3db",
   "metadata": {},
   "source": [
    "### Specify map parameters\n",
    "\n",
    "Now let's tell the map where to center itself! This is done with the \"setCenter\" command applied to the Map object. The first two arguments below are the longitude and latitude, respectively; the third is the zoom level for the map (smaller numbers = bigger maps)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ae32bb3f-1fb0-4a0b-a5bd-2c0b363432c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Map.setCenter(-90.162, 29.8597, 4)   # New Orleans, USA\n",
    "# Map.setCenter(-114.9774, 31.9254, 4) # Mouth of the Colorado River, Mexico\n",
    "# Map.setCenter(-111.1871, 37.0963, 3) # Lake Powell, USA\n",
    "# Map.setCenter(149.412, -35.0789, 4)  # Lake George, Australia\n",
    "Map.setCenter(105.26, 11.2134, 4)     # Mekong River Basin, SouthEast Asia\n",
    "# Map.setCenter(90.6743, 22.7382, 4)   # Meghna River, Bangladesh\n",
    "# Map.setCenter(81.2714, 16.5079, 4)   # Godavari River Basin Irrigation Project, India\n",
    "# Map.setCenter(14.7035, 52.0985, 4)   # River Oder, Germany & Poland\n",
    "# Map.setCenter(-59.1696, -33.8111, 4)  # Buenos Aires, Argentina\\\n",
    "#Map.setCenter(-74.4557, -8.4289, 4)  # Ucayali River, Peru\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d227ff7d-f3c8-4204-b8f2-61a5cdc60832",
   "metadata": {},
   "source": [
    "Now we can set a color palette to use when plotting the data layer. The following are the palettes specified for temperature and precipitation in the GEE description page for ERA5. These are what are called \"Hex\" color codes: don't worry about them a huge amount for now, just know that GEE has a lot of color tables like this that you can look up. \n",
    "\n",
    "Note that you also have to specify a max/min value for each of the variables!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8b6a15e7-061b-4990-833e-a957528cbabb",
   "metadata": {},
   "outputs": [],
   "source": [
    "VIS_TS = {\n",
    "    'min':273,\n",
    "    'max':310,\n",
    "    'palette': ['#000080', '#0000D9', '#4000FF', '#8000FF', '#0080FF', '#00FFFF', '#00FF80',\n",
    "    '#80FF00', '#DAFF00', '#FFFF00', '#FFF500', '#FFDA00', '#FFB000', '#FFA400',\n",
    "    '#FF4F00', '#FF2500', '#FF0A00', '#FF00FF']\n",
    "}\n",
    "\n",
    "VIS_PREC = {\n",
    "    'min':0,\n",
    "    'max':0.1,\n",
    "    'palette': ['#FFFFFF', '#00FFFF', '#0080FF', '#DA00FF', '#FFA400', '#FF0000']\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcc27258-cf68-45ef-a521-51ce2968d526",
   "metadata": {},
   "source": [
    "### Add data to the map\n",
    "\n",
    "Finally, now that we have our precipitation image `prdflt`, we can put it on top of the basemap! This is done using the `.addLayer` function, which we apply to the `Map` object that contains the Google Maps basemap. We also pass it other *arguments*:\n",
    "- the visualization parameters (color ranges, maps) stored in VIS_PREC\n",
    "- the name of the field ('total precipitation')\n",
    "- the *opacity* so that we can also see the basemap underneath (feel free to play with this one to see how things change!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "48bfd420-027d-46f8-b011-15c7402e88b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c8da3429df5f475383fe784c3bb72a67",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(bottom=2220.0, center=[11.2134, 105.26], controls=(WidgetControl(options=['position', 'transparent_bg'], w…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Map.addLayer(prdflt, VIS_PREC,'total precipitation',opacity=0.3)\n",
    "\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c5dba8b-b58b-41d0-9385-64b52df7922e",
   "metadata": {},
   "source": [
    "## Play with the map!\n",
    "\n",
    "Now that everything has been set up, try playing around with things a bit. \n",
    "\n",
    "#### 1) Look at time variations\n",
    "\n",
    "Change the date range in the ee.Filter.date command above, then replot the map layer. How different do things look? \n",
    "\n",
    "#### 2) Change map locations\n",
    "\n",
    "Try zooming in and out on different parts of the map. What types of features do you notice? Are there particular patterns that seem to persist in particular places? (for example: check out the equatorial Pacific!)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43cc9371-ed0c-42cf-86e2-57303e4e8dd9",
   "metadata": {},
   "source": [
    "## Look for historical storm events\n",
    "\n",
    "Test out your map manipulation skills! Try filtering the data in time and space to locate a couple of particularly strong precipitation events:\n",
    "\n",
    "- Montecito, California: January 9, 2018\n",
    "\n",
    "- Rio Grande do Sul, Brazil: March 11-13, 2017\n",
    "\n",
    "- Mumbai, India: July 26, 2005\n",
    "\n",
    "(hint: Google is your friend for finding lats/lons for these places!)\n",
    "\n",
    "Can you see any coherent patterns of precipitation in those locations on those days?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20c2aa59-93b7-4d24-bfa6-a61380963add",
   "metadata": {},
   "source": [
    "## Save images of your maps\n",
    "\n",
    "Finally, save your work! There are a few different options for saving images from GEE; here we'll use the getThumbURL() method. This generates a JPG or PNG image, for which the URL is provided, then you can download the result. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5ce53076-029a-4372-a644-aa31913149c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/1a0117e07448fda243811263823fed20-c54d923cdae299e1253b2d0f747fdb96:getPixels\n"
     ]
    }
   ],
   "source": [
    "url = prdflt.getThumbUrl({\n",
    "    'min': 0, 'max': 0.1, 'dimensions': 512, \n",
    "    'palette': ['#FFFFFF', '#00FFFF', '#0080FF', '#DA00FF', '#FFA400', '#FF0000']})\n",
    "print(url)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62253c3b-5014-43de-977c-4637dcc3835e",
   "metadata": {},
   "source": [
    "Once you've executed the above code, you can save it to whatever filename you like!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a7c7842-4af1-4446-a769-34a2385192d3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}