python-gee.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "61a7ed5d",
   "metadata": {},
   "source": [
    "# Python and Google Earth Engine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50b45947",
   "metadata": {},
   "source": [
    "# To Do!\n",
    "- [ ] Create an exercise\n",
    "- [ ] finish up"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6208c8b",
   "metadata": {},
   "source": [
    "### A note on running this tutorial\n",
    "\n",
    "We recommend using Jupyter Notebook to run the code in this tutorial. The code today makes use of interactive widgets, which will display inline when Jupyter Notebooks. You can run the code in an IDE like spyder, but where the variable `Map` is called, you will have to export the map and look at it in html using the following code: \n",
    "\n",
    "```python\n",
    "Map.to_html(filename='map.html, title='My Map', width='100%', height='880px')\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f79cea7",
   "metadata": {},
   "source": [
    "## Learning goals for today: \n",
    "\n",
    "- Get introduced to object oriented programming and learn to read javascript documentation\n",
    "- Visualize raster and vector data using the google earth engine in python\n",
    "- Use the python api for earth engine for doing geospatial analysis. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dac6b75",
   "metadata": {},
   "source": [
    "Many of you may already be familiar with or have explored Google Earth (https://earth.google.com/web/). Within Google Earth, you can virtually explore the Earth from your computer. The application offers various maps and information, including the familiar Google Maps base map and additional layers such as the age of the sea floor. Moreover, Google Earth (or more accurately, its \"engine\") serves as a powerful cloud-based computational platform for geospatial analysis. It provides a wealth of information and tools that are accessible to everyone (with a Google account). \n",
    "\n",
    "\n",
    "\n",
    "Google Earth Engine (GEE) is the underlying platform behind Google Earth, initially developed and released as a JavaScript API. But what exactly does that mean? Let's break it down into JavaScript and API.\n",
    "\n",
    "### API\n",
    "An API stands for \"Application Programming Interface\" and serves as a standardized language or set of rules that allows different software components to communicate with each other. APIs are commonly used for tasks such as accessing data (e.g., weather data from OpenWeatherMap) or interacting with  [a command line tool: GDAL](https://geoscripting-wur.github.io/PythonRaster/) with python or R, as we covered in the previous tutorials.\n",
    "\n",
    "### JavaScript\n",
    "JavaScript (JS) is a widely used programming language primarily employed for web development. It follows an object-oriented programming (OOP) paradigm, which emphasizes creating reusable code through objects. Objects in JS can encapsulate both data (properties) and functionality (methods) within a single entity. OOP differs somewhat from the programming style we have used thus far, although Python can also be used as an OOP language. Working with the GEE Python API offers a fantastic opportunity to become familiar with object-oriented programming while leveraging the immense power of Google's engine."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec9a1405",
   "metadata": {},
   "source": [
    "# Object-Oriented Programming in Python\n",
    "\n",
    "In this tutorial, we will learn how to interact with GEE. However, before we dive into that, it's essential to understand the principles of Object-Oriented Programming (OOP) and how it is implemented in Python.\n",
    "\n",
    "In Python, objects are created and manipulated using classes. A class serves as a blueprint that defines the structure and behavior of objects. It brings together data (properties) and functions (methods) into a single object.\n",
    "\n",
    "To define a class in Python, we use the `class` keyword, followed by the name of the class. Let's take a look at an example of a simple class called `Person`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e08f760",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Person:\n",
    "    def __init__(self, name, age):\n",
    "        self.name = name  # < this is a property\n",
    "        self.age = age    # < this is also a property\n",
    "    \n",
    "    def greet(self):  # < This \"function\" is a method\n",
    "        print(f\"Hello, my name is {self.name} and I'm {self.age} years old.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6cae57b",
   "metadata": {},
   "source": [
    "In the provided code, the `__init__` method is a special method known as a *constructor*. It is automatically invoked when an object is created from the class. The `self` parameter refers to the instance of the class itself, allowing access to its properties and methods. Whenever a method is defined within a class, we give `self` as the first parameter. \n",
    "\n",
    "The `Person` class has two properties, `name` and `age`, as well as one method, `greet`. The `greet` method prints a greeting message that includes the person's name and age.\n",
    "\n",
    "To create an instance of the `Person` class, you simply call the class as if it were a function and assign the result to a variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4b5e90e",
   "metadata": {},
   "outputs": [],
   "source": [
    "person1 = Person(\"Alice\", 25)\n",
    "person2 = Person(\"Bob\", 30)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8a320a6",
   "metadata": {},
   "source": [
    "We have created two objects, `person1` and `person2`, which are instances of the `Person` class. Now we can access the properties and call the methods of these objects:\n",
    "\n",
    "```python\n",
    "print(person1.name)  # Output: Alice\n",
    "print(person2.age)   # Output: 30\n",
    "person1.greet()      # Output: Hello, my name is Alice and I'm 25 years old.\n",
    "person2.greet()      # Output: Hello, my name is Bob and I'm 30 years old.\n",
    "```\n",
    "\n",
    "This example is straightforward, but keep in mind that classes can become more complex. For instance, take a look at the creation of a GeoPandas `GeoSeries` [here](https://github.com/geopandas/geopandas/blob/main/geopandas/geoseries.py). It may seem complicated and hard to read, but it is well-documented, so try to understand some of its functionality. The `to_json` method, for example, provides straightforward functionality (though it may still be challenging to read due to its complexity).\n",
    "\n",
    "You may have noticed that class definitions differ from what we have learned. In the GeoPandas example, the class is defined as follows:\n",
    "\n",
    "```python \n",
    "class GeoSeries(GeoPandasBase, Series):\n",
    "```\n",
    "\n",
    "However, we previously learned to define a class like this:\n",
    "\n",
    "```python \n",
    "class GeoSeries:\n",
    "```\n",
    "\n",
    "The difference lies in the code within the parentheses, which represents *inheritance*. The class will inherit all the functionality from the classes specified as arguments, in this case, `GeoPandasBase` and `Series`. The `Series` object refers to the [Pandas.Series](https://github.com/pandas-dev/pandas/blob/main/pandas/core/series.py#L243C15-L243C15), which contains thousands of lines of code with various functionalities. Therefore, the GeoPandas `GeoSeries` contains all the functionalities implemented in Pandas, as well as those from `GeoPandasBase` and `GeoPandas.Series` itself.\n",
    "\n",
    "Phew, that's a lot of complicated code, and it may seem overwhelming. You don't need to understand every detail. The important thing is to grasp the value of classes, objects, and inheritance. How does this work in a simple example?\n",
    "\n",
    "So, we've learned that with inheritance, classes can acquire all the functionality from other classes and build upon that. Let's create a new object called `Student`, which will inherit all the properties and methods from the `Person` class we created earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8916711",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Student(Person):\n",
    "    def __init__(self, name, age, student_id):\n",
    "        super().__init__(name, age)\n",
    "        self.student_id = student_id\n",
    "    \n",
    "    def study(self):\n",
    "        print(f\"{self.name} is studying.\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "693378b2",
   "metadata": {},
   "source": [
    "In the provided code, the `Student` class inherits from the `Person` class, which we will refer to as the superclass. By doing so, it extends the functionality of the superclass by adding a new property (`student_id`) and a new method (`study`). The `super()` function is used to call the superclass's `__init__` method, allowing the subclass to initialize the inherited properties.\n",
    "\n",
    "As a result, the `Student` class contains both the methods and properties inherited from the `Person` class, as well as the additional ones defined within the `Student` class:\n",
    "\n",
    "```python\n",
    "student = Student(\"Eve\", 22, \"123456\")\n",
    "print(student.name)         # Output: Eve\n",
    "print(student.student_id)   # Output: 123456\n",
    "student.greet()             # Output: Hello, my name is Eve and I'm 22 years old.\n",
    "student.study()             # Output: Eve is studying.\n",
    "```\n",
    "\n",
    "In this example, `student` is an instance of the `Student` class. It can access the inherited properties from the `Person` class, such as `name`, as well as the newly added property `student_id`. Similarly, it can invoke both the inherited method `greet` and the additional method `study`, which are specific to the `Student` class."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f39d27f",
   "metadata": {},
   "source": [
    "# Google Earth Engine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0003c64",
   "metadata": {},
   "source": [
    "Google Earth Engine (GEE) is a comprehensive platform that offers a vast collection of geospatial data along with powerful computational capabilities. It enables users to analyze and visualize geospatial data on a global scale, making it a valuable tool for various research domains such as environmental monitoring and resource management.\n",
    "\n",
    "One of the prominent features of GEE is its extensive collection of raster data sources. You can explore the catalog of available datasets [here](https://developers.google.com/earth-engine/datasets/catalog). Additionally, GEE provides access to raster time series and non-raster data as well. For specific use cases and examples, you can refer to the article by Pérez-Cutillas and colleagues published this year [here](https://www.sciencedirect.com/science/article/pii/S2352938522002154). They conducted a systematic review of research conducted using GEE. From their article:\n",
    "\n",
    "\n",
    "> The results of the meta-analysis following the systematic review showed that: (i) the Landsat 8 was the most widely-used satellite (25%); (i) the non-parametric classification methods, mainly Random Forest, were the most recurrent algorithms (31%); and (iii) the water resources assessment and prediction were the most common methodological applications (22%).\n",
    "\n",
    "-- Pérez-Cutillas et al., 2023"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31a6a6af",
   "metadata": {},
   "source": [
    "## Nighttime Light Emission\n",
    "\n",
    "In this tutorial, we will conduct a relatively straightforward analysis to demonstrate commonly used tools in Google Earth Engine (GEE). Our main focus is to show you how to work with GEE in Python, locate tools and documentation, and read the JavaScript (JS) documentation.\n",
    "\n",
    "The analysis we will perform involves comparing the amount of light emitted at night in the Netherlands over time. However, you can choose a different area for your analysis, as we will be using a global dataset from 2012 to the present.\n",
    "\n",
    "To set up the required environment, you will need the following dependencies. Save the following YAML code in a new file named `env.yaml` and create and activate the environment using mamba \n",
    "\n",
    "```yaml\n",
    "name: pythonGee\n",
    "channels:\n",
    "  - conda-forge\n",
    "dependencies:\n",
    "  - python=3.10\n",
    "  - jupyter\n",
    "  - earthengine-api\n",
    "  - geemap\n",
    "  - pygadm\n",
    "  - geopandas\n",
    "  - geojson\n",
    "```\n",
    "\n",
    "To test the environment, try importing the following packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d11200fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import geemap \n",
    "import pygadm\n",
    "import geopandas as gpd\n",
    "import geojson"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5e84256",
   "metadata": {},
   "source": [
    "To work with Google Earth Engine (GEE), you will need a Google account for authorization. The authentication process can be initiated by calling the `Authenticate` function from the `ee` module. Follow these steps:\n",
    "\n",
    "1. After running the `Authenticate` function, a URL will be printed in the output.\n",
    "1. Click on the URL, which will redirect you to a Google account login page.\n",
    "1. Log in to your Google account to create an authorization token.\n",
    "1. When prompted with a warning stating that Google has not verified the app, click \"Continue\". This warning is shown because the code we are about to run is not verified by Google, ensuring your account's security.\n",
    "1. On the authorization page, select all the checkboxes. This grants permission to access GEE data and manage data and permissions in the cloud storage where the analysis will be performed.\n",
    "1. Finally, you will see an authorization code. Copy the code and paste it into the text box that appears after running the `Authenticate` function. Press Enter to complete the authentication process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1629c477",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ee.Authenticate()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1121b4ce",
   "metadata": {},
   "source": [
    "You now have (hopefully) succesfully authenticated. Next we have to initialize using the `Initialize()` method: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e877a522",
   "metadata": {},
   "outputs": [],
   "source": [
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b0d88e1",
   "metadata": {},
   "source": [
    "If you do not encounter any warnings or errors, you have successfully authenticated and are now ready to work with the Earth Engine library in Python! Congratulations! \n",
    "\n",
    "The next time you work in the same environment, you won't need to authenticate again. Simply run the `Initialize` function, and you'll be all set to use the Earth Engine library. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67531ce7",
   "metadata": {},
   "source": [
    "## Some notes on documentation \n",
    "\n",
    "In the documentation we will refer to today, you will notice that there are examples available for both JavaScript and Python (Google Colab). This is because Google is working on implementing all the JavaScript functionality in Python. Not all the documentation or examples are at this time available in Python. Luckily, although the languages are quite different, the JavaScript and Python implementation for this package are very similair. It is a good skill to be able to read some JavaScript. Additionally, we are not sure what documentation is available for python at the time of following this tutorial. Therefor we will refer to the JavaScript documentation, but we will highlight any important differences and provide guidance on using this code in Python.\n",
    "\n",
    "Some common differences are the following. To initialize a variable, in Python you can simply type: \n",
    "\n",
    "```python\n",
    "variablename = 'value'\n",
    "```\n",
    "\n",
    "In Javascript we can initialize a variable in various ways. The difference is where this variable in the code will be known: \n",
    "\n",
    "```Javascript\n",
    "var variablename = 'value'\n",
    "let variablename = 'value'\n",
    "const variablename = 'value'\n",
    "```\n",
    "\n",
    "Another often occuring difference in the documentation is how variables are passed to a function. In python we can do that as we are used to: \n",
    "\n",
    "```python\n",
    "multiStats = img.reduceRegions(\n",
    "  collection=regionCol,\n",
    "  reducer=reducer,\n",
    "  scale=30,\n",
    "  crs='EPSG:3310',\n",
    ")\n",
    "```\n",
    "\n",
    "In Javascript we have to pass an object, similar to a dictionary in python:\n",
    "\n",
    "```javascript\n",
    "var multiStats = img.reduceRegions({\n",
    "  collection: regionCol,\n",
    "  reducer: reducer,\n",
    "  scale: 30,\n",
    "  crs: 'EPSG:3310',\n",
    "});\n",
    "```\n",
    "\n",
    "You do not need to understand this code itself for now, just note the differences in syntax. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a83b528b",
   "metadata": {},
   "source": [
    "## Accessing Night Light Data\n",
    "\n",
    "To analyze the development of light emission over the past 10 years, we will compare the amount of emitted light per province on a monthly basis.\n",
    "\n",
    "For this analysis, we will need data on emitted light. We will use the VIIRS Nighttime data, which you can learn more about [here](https://developers.google.com/earth-engine/datasets/catalog/NOAA_VIIRS_DNB_MONTHLY_V1_VCMCFG). This dataset provides monthly average radiance composite images using nighttime data from the Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB). You can also explore the data in [this Esri viewer](https://www.arcgis.com/apps/mapviewer/index.html?layers=edabcbb5407547f5bc883018eb6e7986). To work with this data, we can create an `ImageCollection` object. You can refer to the [documentation](https://developers.google.com/earth-engine/apidocs/ee-imagecollection) for more details.\n",
    "\n",
    "\n",
    "\n",
    "To create an `ImageCollection`, you can use the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3099281e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the night time light emision data.\n",
    "viirs_image_collection = ee.ImageCollection(\"NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c2a264e",
   "metadata": {},
   "source": [
    "Note that this is the exact same code as in the JS documentation, but adapted for this tutorial. In Python, to retrieve information about an object in the form of a dictionary, we need can call the `getInfo()` method.\n",
    "\n",
    "```JavaScript\n",
    "print('Image collection from a string', ee.ImageCollection(\"NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG\"));\n",
    "```\n",
    "\n",
    "In Python, we use the following command:\n",
    "\n",
    "```python\n",
    "print(ee.ImageCollection(\"NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG\").getInfo()['properties'])\n",
    "```\n",
    "\n",
    "Here, we call the `getInfo` method to convert the object information to a dictionary that we can extract useful information from. However, when working in Jupyter notebooks, we can simply call the object itself to display it inline:\n",
    "\n",
    "```python\n",
    "ee.ImageCollection(\"NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG\")\n",
    "```\n",
    "\n",
    "This will display the ImageCollection object in Jupyter. \n",
    "\n",
    "Inspecting the ImageCollection object, we can see that it has several properties. It includes `type`, `id`, `version`, an empty list for `bands`, 23 `properties`, and 132 `features`. These properties are similar to what you would expect in a [GeoJSON](https://geojson.org/) representation of a FeatureCollection.\n",
    "\n",
    "The `type`, `id`, and `version` provide information about the data and the object type. At this point, the object has an empty list for bands, indicating that there are no bands defined yet. It has 23 properties, which describe information about the collection as a whole (e.g., it is a monthly collection). The collection contains features, which are individual images. Each feature (image) has its own `type`, `id`, `version`, `bands`, and `properties`. There are no features within these images themselves, as are no collections.\n",
    "\n",
    "\n",
    "## Earth Engine as a viewer\n",
    "Now, let's explore how we can use Google Earth Engine as a viewer. As mentioned earlier, Google Earth Engine is not only a data collection but also a platform with tools and a viewer. To interact with the viewer, we will utilize a package called `geemap`. Geemap also provides useful tools, such as data conversion to and from GEE native objects. We'll cover this in more detail later.\n",
    "\n",
    "To create a map, we will follow a few steps outlined in the [geemap getting started page](https://geemap.org/get-started/#add-earth-engine-data), modifying it for our use case:\n",
    "\n",
    "1. Initialize a map with a starting position and zoom level. This map will serve as our base, and we will add layers to it.\n",
    "1. Select one image from our image collection.\n",
    "1. Add the selected image to the map.\n",
    "1. Add a boundary of the Netherlands to the same map.\n",
    "\n",
    "\n",
    "#### Initialize the map \n",
    "\n",
    "By following these steps, we will create a map visualization of the selected image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ab0a48d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set our initial zoom position for the netherlands\n",
    "center_lon = 51.962589\n",
    "center_lat = 5.669627\n",
    "zoomlevel = 8\n",
    "\n",
    "# The following line will initialize the map\n",
    "Map = geemap.Map(center=[center_lon, center_lat], zoom=zoomlevel)\n",
    "\n",
    "# The following line will actually call and show the map:\n",
    "Map\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "987b2d5c",
   "metadata": {},
   "source": [
    "The map will stay and remain in this window unless we instruct it otherwise. It is an interactive map with 1 standard basemap at the moment. To start over, you can call `Map.clear()`, which will remove all layers and other objects from the map, requiring you to initialize it again. Sometimes, when working with an interactive widget like this map, Jupyter might become unresponsive. If you no longer receive any response (you can test this by running `1 + 1` in a new code block), you might need to restart the kernel and begin again.\n",
    "\n",
    "\n",
    "#### Select 1 image and add it to the map \n",
    "\n",
    "Step 1 is complete! The next step is to select one image from the image collection and add it to the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5b6c465",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial date of interest (inclusive).\n",
    "from_date = '2017-01-01'\n",
    "# Final date of interest (exclusive).\n",
    "to_date = '2023-01-01'\n",
    "\n",
    "image = viirs_image_collection\\\n",
    "            .filterDate(from_date, to_date)\\\n",
    "            .sort(\"system:time_start\", True)\\\n",
    "            .first()\n",
    "\n",
    "visualization_params = {\n",
    "    'min': 0,\n",
    "    'max': 10,\n",
    "    'opacity': 0.5,\n",
    "    'bands': ['avg_rad']\n",
    "}\n",
    "\n",
    "Map.addLayer(ee_object = image, \n",
    "             vis_params = visualization_params,\n",
    "             name = 'My First Image')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b992f156",
   "metadata": {},
   "source": [
    "In the code snippet provided, several steps are performed to add an image to the map. In the first line the following happened:\n",
    "\n",
    "1. Filtering Images: The `viirs_image_collection` is filtered to select only the images between two specified dates using the `filterDate(from_date, to_date)` method.\n",
    "\n",
    "1. Sorting the Collection: The filtered image collection is sorted based on the `\"system:time_start\"` property in ascending order using the `sort(\"system:time_start\", ascending)` method. The resulting collection is then assigned to the variable `image`.\n",
    "\n",
    "1. Selecting an image: We selected the first image out of the collection using the `first()` method. \n",
    "\n",
    "We then define visualization parameters and add the image to the map: \n",
    "\n",
    "1. Visualization Parameters: Visualization parameters are defined to specify how the image should be displayed on the map. This includes the desired range of values, opacity, and the specific band to visualize (`avg_rad` in this case).\n",
    "\n",
    "1. Adding the Image to the Map: The `image` is added to the map using the `AddLayer` method. This method requires the Earth Engine object (`image`), the visualization parameters, and an optional name for the layer.\n",
    "\n",
    "By following these steps, the code successfully adds the selected image to the map, great!\n",
    "\n",
    "#### Add the boundaries of the Netherlands to the map\n",
    "\n",
    "The next code will add the boundaries of the Netherlands to the map. To get the boundary data we use the `pygadm` package. `pygadm` is a dataset project from the University of Berkeley that provides administrative boundary data on multiple scales for various regions worldwide. Be careful! The returned data is not always completely correct, but it suits our purposes for now. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a99d4db",
   "metadata": {},
   "outputs": [],
   "source": [
    "netherlands = pygadm.get_items(admin='NLD', content_level=1)\n",
    "ee_netherlands = geemap.geopandas_to_ee(netherlands)\n",
    "\n",
    "Map.addLayer(ee_object = ee_netherlands, \n",
    "             name = 'Netherlands Boundaries')\n",
    "Map.addLayerControl()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d69b0bc0",
   "metadata": {},
   "source": [
    "Using the `get_items` function, we can obtain the boundaries at content level 1, which represent the provinces in the Netherlands. This data is retrieved in the `GeoPandas GeoDataFrame` format, which cannot be added to the map directly. However, we can utilize a convenient conversion tool provided by the `geemap` package to convert it into an Earth Engine feature. The resulting object can then be easily placed on the map.\n",
    "\n",
    "For additional examples and tools, we recommend checking out the [examples](https://github.com/gee-community/geemap/tree/master/examples/notebooks) provided by `geemap` themselves! It's worth noting that you can also [export a map to HTML](https://github.com/gee-community/geemap/blob/master/examples/notebooks/21_export_map_to_html_png.ipynb), which might be useful for your project.\n",
    "\n",
    "Take a look at these examples to further enhance your project and explore the capabilities offered by `geemap`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "778731c7",
   "metadata": {},
   "source": [
    "Great, we have now made a simple map in google earth engine. We added an image for the whole world and we imported external data and also visualized that. And it works very quick!  \n",
    "\n",
    "## Analysis in GEE \n",
    "\n",
    "For the analysis we want to clip the imagery to the boundaries. After that we want to take the average of all the pixels per province. This will result in a single value per province per month. Schematically the steps will look like this:\n",
    "\n",
    "1. Select data: from the ImageCollection, select the band we want to do the analysis on and filter the collection to the time range we want to do the analysis on.\n",
    "1. Single image: Collapse the image collection into a single image with many bands\n",
    "1. Clip: Clip the images on the boundaries\n",
    "1. Calculate the statistics per province\n",
    "\n",
    "#### Selecting the data and collapsing it\n",
    "\n",
    "The first step can be done in one line. Several times we call a method on an image collection, which will result in an image collection. On this collection we will call another method and so fort. Look at the following line of code and try to understand what is happening before you read the explanation. Use the documentation!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6259078",
   "metadata": {},
   "outputs": [],
   "source": [
    "viirs_image = viirs_image_collection.select('avg_rad') \\\n",
    "                .filterDate(from_date, to_date) \\\n",
    "                .toBands() "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "183216ac",
   "metadata": {},
   "source": [
    "4 methods are called in one line. In this code we go from an image collection over all available time steps to a single image of the time we are interested in, with one band representing an image per month. Why we collapse it into this single image will become clear in the next step. The separate methods are described below.\n",
    "1. First, we select the band we are interested in, the `avg_rad` band, containing the avarage nighttime radiation per month. For this we use the [`select`](https://developers.google.com/earth-engine/apidocs/ee-imagecollection-select) method. \n",
    "1. We then select the images from the dates we are interested in using the [`filterDate`](https://developers.google.com/earth-engine/apidocs/ee-imagecollection-filterdate) method. \n",
    "1. Finally we call the [`toBands`](https://developers.google.com/earth-engine/apidocs/ee-imagecollection-tobands) method. Note that this method results in an *`Image`*  instead of an *`ImageCollection`*. \n",
    "\n",
    "When you run the code above note that the code is run instantly, no matter the processing power of your laptop or virtual machine, while actually it is quite a large computation. How can this be? Actually the code and the calculation is not done on your computer, in fact there is no computing done at all yet for this code. This code creates an API Call, which will be sent to the Earth Engine cloud platform once it needs to be called. For example when you now call the `viirs_image` variable, it takes longer. The API call is sent and the results are returned to your environment.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8fa6837",
   "metadata": {},
   "outputs": [],
   "source": [
    "viirs_image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5230593",
   "metadata": {},
   "source": [
    "#### Clipping the data\n",
    "\n",
    "Next, we will clip the image to the desired area, the boundary of the Netherlands. This can be done by using the [`clip`](https://developers.google.com/earth-engine/apidocs/ee-image-clip) method. We could have also added this to the already long line of code above. Functionally, this makes no difference. The decision on how to do this is personal; every programmer might do it differently. The most important aspect of this choice is to make it clear for yourself and others what each line is doing and to define clear variable names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc85057d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to clip an image to the specified geometry\n",
    "viirs_image_clipped = viirs_image.clip(ee_netherlands)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a846d77",
   "metadata": {},
   "source": [
    "`viirs_image_clipped` now contains an image with 36 bands, verify this in your own environment. Do you understand the band names? \n",
    "\n",
    "We now want to calculate the average per province. We do this using a *reducer*, which we apply over an area. A reducer is an algorithm which will be applied over a set of values. The set of values in our case will be the pixel values within a province. How this works in earth engine, is that we use a method called [`reduceRegions`](https://developers.google.com/earth-engine/apidocs/ee-image-reduceregions). According to the documentation this method needs as input a collection and a reducer. As input we will give the polygons of the provinces in the Netherlands: `ee_netherlands`. As reducer we have options. To see what the options are, search for `reducer` in the filter in the documentation. Since we want the average over the collection, we will use the `ee.Reducer.mean()` function, which is not a reducer itself but it *returns* a reducer. That means we will need to actually call the function in the input! We also pass give the optional parameter scale, and we set it to 300 meters. Check the documentation to understand what this parameter does. \n",
    "\n",
    "`viirs_image_clipped` now contains an image with 36 bands. Please verify this in your own environment. Do you understand the band names?\n",
    "\n",
    "#### Calculating statistics\n",
    "\n",
    "Now, we want to calculate the average radiance per province. We can achieve this using a *reducer*, which we apply over an area. A reducer is an algorithm that is applied over a set of values. In our case, the set of values will be the pixel values within a province. We will use a method called [`reduceRegions`](https://developers.google.com/earth-engine/apidocs/ee-image-reduceregions) to gather this set of values. According to the documentation, this method requires a collection and a reducer as inputs. We will provide the polygons of the provinces in the Netherlands, referred to as `ee_netherlands`, as the input collection. As for the reducer, there are various options available. To explore the available options, you can search for `reducer` using the filter in the documentation.\n",
    "\n",
    "Since we want to take the average over the collection, we will utilize the `ee.Reducer.mean()` function. It is important to note that this function itself is not a reducer, but rather it *returns* a reducer. Therefore, we need to actually call the function as part of the input. Additionally, we provide the optional parameter `scale` with a value of 285 meters. \n",
    "\n",
    "Why do we chooe 285 meters as scale? Find out using the documentation of the imagery and the documentation of `reduceRegions`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "189d2b98",
   "metadata": {},
   "outputs": [],
   "source": [
    "province_statistics = viirs_image_clipped.reduceRegions(\n",
    "    collection = ee_netherlands,\n",
    "    reducer = ee.Reducer.mean(),\n",
    "    scale = 285\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8007498d",
   "metadata": {},
   "source": [
    "## Doing external analysis and visualization\n",
    "\n",
    "All right! We now have a FeatureCollection with the average radiation per province, per month over the given time period. To perform analysis and visualization, our next step is to calculate the largest difference within the time period, classify this difference, and display the results on a map. While this is all possible within Earth Engine, since you are familiar with working with GeoPandas, it makes sense use what we know. Exporting the output to a GeoPandas DataFrame and continuing from there is not a problem at all, so let's do that.\n",
    "\n",
    "The steps we will follow are as follows:\n",
    "1. Export the statistics FeatureCollection to a GeoPandas DataFrame.\n",
    "1. Calculate the difference in radiation per province.\n",
    "1. Classify this difference into three classes: small difference, medium difference, and large difference.\n",
    "1. Convert the resulting DataFrame back into an Earth Engine object.\n",
    "1. Visualize the results.\n",
    "\n",
    "\n",
    "#### Exporting\n",
    "`Geemap` provides useful tools for this task. The most convenient format for us is a [GeoPandas](https://geemap.org/common/?h=ee_to_geotiff#geemap.common.ee_to_geopandas) object, but there are many other options available as well. Please refer to the documentation for additional options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad8be9a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = geemap.ee_to_geopandas(province_statistics)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6645ba79",
   "metadata": {},
   "source": [
    "#### Calculating differences\n",
    "\n",
    "From here calculating the difference is straight forward, since we know (geo)pandas well. We now also have all the functionality that comes with GeoPandas, such as very easily writing the results to a geojson!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "054b67ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['diff'] = df['20170101_avg_rad'] - df['20230101_avg_rad']\n",
    "df.to_file('test.geojson', drive='GeoJSON')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a2a2bee",
   "metadata": {},
   "source": [
    "#### Classifying and visualizing trends\n",
    "\n",
    "Now to see how the radiation behaved over time we have to do a little more work, since the dataframe in this format is not very useful for this. Read and understand the code below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c76c9db8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2cbbb8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean the dataframe a little\n",
    "clean_df = df.drop(['GID_0', 'GID_1', 'ISO_1', 'HASC_1', 'NAME_0', 'NL_NAME_1', 'VARNAME_1', 'TYPE_1','CC_1', 'ENGTYPE_1'], axis=1)\n",
    "\n",
    "# Reformat the dataframe so the rows depict measurements per timestamp per province, instead of an attribute per timestamp\n",
    "melted_df = pd.melt(clean_df, id_vars=['NAME_1', 'geometry', 'diff'] , var_name='date', value_name='value')\n",
    "\n",
    "# Convert the time string to a datetime format\n",
    "melted_df['short_date'] = melted_df['date'].str[:8].str.strip()\n",
    "melted_df['datetime'] = pd.to_datetime(melted_df['short_date'], format=\"%Y%m%d\")\n",
    "\n",
    "\n",
    "# Initialize the plot\n",
    "fig, ax = plt.subplots(figsize=(9, 6))\n",
    "\n",
    "# Create an empty list to save all lines in for the legend \n",
    "lines = []\n",
    "\n",
    "# For each province, add a line with a different color \n",
    "for prov in melted_df['NAME_1'].unique():    \n",
    "    \n",
    "    # Create a dataframe for the province of interest\n",
    "    d = melted_df[melted_df['NAME_1'] == prov]\n",
    "    \n",
    "    # Generate a random color\n",
    "    random_rgb_values = list(np.random.random(3))\n",
    "    random_color = matplotlib.colors.rgb2hex(random_rgb_values, keep_alpha=False)\n",
    "    \n",
    "    # Add a line for each province\n",
    "    lines += ax.plot( d[\"datetime\"].values, d[\"value\"].values, color= random_color, lw=1, zorder=10, label=prov)\n",
    "    \n",
    "    \n",
    "# Add the legend    \n",
    "ax.legend(handles = lines)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39d202ec",
   "metadata": {},
   "source": [
    "We can see that there are some differences per province, but that there does not seem to a nationwide trend. We can also see that in the summer there is no data for the Netherlands, or in fact for entire northern europe. This has to do with the nights being so long here in the summer. Per province there might be a larger difference. To find out, we can add a trendline to the graph, but this is out of scope for this tutorial. \n",
    "\n",
    "#### Classify differences\n",
    "For the visualization we now want to categorize the differences in 3 classes, a little difference, medium difference and a large difference. To do this we will use a numpy function called `linspace` and a pandas function called `cut`. This results in an array with datatype `CategorizalClasses`, which earth engine does not understand so we convert it to a string. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7297867",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create 3 classes\n",
    "breaks = np.linspace(df['diff'].min() - 0.0001, df['diff'].max(), 4)\n",
    "\n",
    "# Assign values to classes\n",
    "classes = ['Low', 'Mid', 'High']\n",
    "df['class'] = pd.cut(df['diff'], bins=breaks, labels=classes, right=True)\n",
    "\n",
    "# Convert the datatype to strings so earth engine understands\n",
    "df['class'] = df['class'].astype('str')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05794d8c",
   "metadata": {},
   "source": [
    "#### Show results on a map\n",
    "\n",
    "So the next step is to show the differences on an earth engine map. To do this, we have to convert the geopandas dataframe back to a earth engine object. Before we do this we have to tell geopandas what coordinate refrence system the geometries are in. Strangely the coordinate system is not exported from earth engine, but it does need it to be assigned to be imported... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1be4ea1",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_crs(4326)\n",
    "new_ee_statistics = geemap.geopandas_to_ee(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e62bb9bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Temp\n",
    "new_ee_statistics = geemap.geojson_to_ee('ee_stats.json')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a88fbae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you run the entire tutorial in one go, the map will show at the same location as it was first initialized.\n",
    "# To show it hear, clear the map and create a new one. \n",
    "# If the map was not initizalized before, this line will throw an error. \n",
    "Map.clear()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7099e1d1",
   "metadata": {},
   "source": [
    "#### Style the map\n",
    "\n",
    "Now we have the polygons back in the right format with the classes. To visualize we define a palette, which tells earth engine which color we will give to each class. For the legend we can use this dictionary directly. Of course we also want to add the latest night time imagery behind it and we want to add a layer control. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9249bed",
   "metadata": {},
   "outputs": [],
   "source": [
    "Map = geemap.Map(zoomlevel=8)\n",
    "\n",
    "palette = {\n",
    "    'High': 'B81D13', # red\n",
    "    'Mid': 'EFB700',  # yellow\n",
    "    'Low': '008450'   # green\n",
    "}\n",
    "\n",
    "\n",
    "visualization_params = {\n",
    "    'min': 0,\n",
    "    'max': 10,\n",
    "    'opacity': 0.7,\n",
    "    'bands': ['20230101_avg_rad']\n",
    "}\n",
    "\n",
    "Map.addLayer(ee_object = viirs_image_clipped, \n",
    "             vis_params = visualization_params,\n",
    "             name = 'Night Time')\n",
    "\n",
    "Map.add_styled_vector(new_ee_statistics, column=\"class\", palette=palette, layer_name=\"Styled vector\")\n",
    "Map.center_object(new_ee_statistics)\n",
    "\n",
    "Map.add_legend(legend_title=\"Legend\", legend_dict=palette, position = 'topright')\n",
    "\n",
    "Map.addLayerControl()\n",
    "\n",
    "Mapg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2275954",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "That marks the end of todays tutorial. Today you learned about many different things. We began with an introduction about object oriented programming. We have learned what objects are with methods and properties and how these can be inherited from other classes. We did this so we understand JavaScript a little better and we can now understand the documentation of the Earth Engine API better. \n",
    "\n",
    "We then used the Python implementation of this API to create maps in Earth Engine and we did some analysis about light radiation in the Netherlands. We exported and imported earth engine objects to other formats we can easier work with, so we can use the tools we know in combination with the power offered by Google Earth Engine. "
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}