For those responsible for defining adaptation strategies, it is crucial to have clear information about climate change impacts, consider uncertainties, and understand the magnitude of potential effects. The Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) provides a framework for consistently projecting the impacts of climate change across different sectors and spatial scales. Through the collaboration of an international network of climate-impact modelers, ISIMIP offers a detailed and consistent picture of the world under different climate change scenarios.

General Objectives

This study aims to test the applicability and replicability of the workflow developed by Müller & Döll (2024) to assess the impact of climate change on a specific variable available in the ISIMIP repository: groundwater recharge. Our focus is on an area of interest (AOI) around Lyon, France.

Data Utilized

The study utilized an extensive dataset comprising 212 NetCDF files, amounting to 60 GB of data. This dataset includes:

• Historical Data: 54 NetCDF files (22 GB) covering the period from 1850 to 2014.
• Future Projections: 158 NetCDF files (38 GB) covering the period from 2015 to 2100.
• Global Climate Models (GCMs): Six models were considered — ipsl-cm6a-lr, mri-esm2–0, mpi-esm1–2-hr, ukesm1–0-ll, and gfdl-esm4.
• Global Hydrological Models (GHMs): Four models were included — cwatm, watergap2–2e, h08, and miroc-integ-land.
• Shared Socioeconomic Pathways (SSPs): Three scenarios were considered — SSP 1, SSP 3, and SSP 5.
• ISIMIP 3b Protocol: The most recent and representative protocol was used for this study.

Method

The method of analysis is based on the work described by Müller & Döll (2024):

• Determine a reference distribution of the annual mean groundwater recharge for each model combination over a historical period from 1985 to 2015.
• Define horizons of interest and 30-year distributions (aroun 2050, 2070 and 2084) and calculate the annual mean groundwater recharge for each model combination, then calculate the relative change compared to its reference mean for different SSPs.
• Analyze the interannual variability to determine the relative changes for dry, normal, and wet years. The groundwater recharge values are sorted by magnitude, and relative differences have been calculated between values of the same rank in the future period and the reference period. This sorting corresponds to the exceedance probability, which is indicated as a percentage.

In the following, the results are described by means of percentile boxes (Figure 3) to describe the distribution of results given by the multi-models ensemble.

Results

The analysis based on Figure 4, revealed significant insights into the future of groundwater recharge around Lyon:

• Short-Term Projections (by 2050): Approximately 90% of the models indicate a decrease in groundwater recharge. This trend is observed in both optimistic and pessimistic scenarios.
• Long-Term Projections (by 2084): In the pessimistic scenario (SSP 5), 80% of the models project a decrease in mean groundwater recharge by 10% to 40%.

Additionnaly, the future interannual fluctuations have been analysed for the horizon 2084 and considering the pessimistic scenario (SSP 5)(Figure 5):

• Wet Years (exceedence between 0 to 20%): 80% of the models indicate a decrease in recharge by 10% to 30%.
• Normal Years (exceedence between 40 to 60%): 80% of the models indicate a decrease in recharge by 10% to 40%.
• Dry Years (exceedence between 80 to 100%): 80% of the models indicate a decrease in recharge by 0% to 55%.

Conclusion

The study successfully developed a reproducible and scalable workflow using ISIMIP repository data and the Müller & Döll (2024) methodology to determine the impact of climate change on a complex environmental variable — groundwater recharge — around a location of interest. This workflow is not limited to Lyon or groundwater recharge; it can be applied to any location worldwide and to various other hydrologic or environmental variables, such as agriculture, biodiversity, and marine ecosystems.

The power of this workflow lies in its comprehensive approach, utilizing all available models of the ISIMIP platform, scenarios, and combinations of GCMs, GHMs, and SSPs. By leveraging this extensive dataset, we provide a holistic view of the potential consequences of climate change on variables of interest. This approach offers valuable insights for policymakers and stakeholders, allowing them to make informed decisions and develop robust adaptation strategies that consider a wide range of potential outcomes — from optimistic to pessimistic scenarios.

Acknowledgements

Thanks to Laura Müller, Hannes Müller Schmied and Petra Döll for scientific advices.

Complementary notes — contact

This article and the code has been elaborated with the support of Ageoce. Ageoce builds digital solutions and services focusing on geodata and geosciences.

Please feel free to contact us anytime.

In France, about 10 millions of houses are exposed to the shrinkage and swelling of clay soils [1] and between 1989 and 2021, the cost of damages related to this phenomenon represented more than 15 billion euros [2]. This phenomenon mostly occurs after severe draught events, that can be identified by low values of the soil wetness index (SWI).

The SWI represents the state of the soil water reserve (2 meters deep) in relation to the available water at the field capacity (water available for plant feeding). For more details, see the mathematical definition in appendix 1. This index is a key component of the national compensation progam: when damages are caused by the shrinkage and swelling of clay soils, the official recognition of “natural disaster” relies on the SWI [3]. If the return period of the SWI value is higher than 25 years, then municipalities can be recognized and supported by the national natural disasters compensation program.

In the current context of climate change, the following questions arise:

• How did the SWI evolve in France since the beginning of its calculation at a national scale in 1969?
• Is there eventually some locations more impacted than others in terms of soil wetness evolution?

Dataset availability and description

We explore the evolution of the uniform SWI (see appendix 1) in France between 1970 to 2021. Uniform SWI values are provided by Meteo-France on a monthly basis since 1969 with a spatial resolution of 8 km. Each monthly value integrates the moisture status of the current month and of the two previous months. The dataset is available on the website of Meteo-France.

Google Earth Engine Analysis

We first ingested this dataset in Google Earth Engine as an Image Collection. This Earth Engine asset has been made publicly available.

Then, the dataset has been reduced and resampled regarding the four seasons of interest to eventually differentiate spring, summer, autumn, or winter soil wetness evolutions over the period 1970–2021.

The uniform SWI image collection ingested in Earth Engine is composed by 624 images (1 image by month from January 1970 to December 2021). For each year, and for each season, a composite image has been create with the mean SWI values of the season.

Then, a linear fit reduction has been applied to each seasonal collection in order to get the general trend. This process computes the least squares estimate of a linear function of the seasonal uniform SWI variable.

Finally, the resulting national trend image is reduced at the scale of French departments. It allows for estimating the averaged SWI gain/loss over the period of interest.

Results

The result of the analysis is illustrated on the figure below. On this figure, the absolute SWI gain/loss between 1970–2021 is mapped. The quantitative table associated with these results is provided in appendix 2, where the averaged SWI gain/loss by department and by season is given.

These results can be summarized as follow:

• The 5 most impacted departments in terms of SWI loss are Gers, Cantal, Haute-Garonne, Lozère and Loire.
• The spring is the season where the most severe and generalized SWI loss is observed. The most severe spring SWI loss is observed in Ardèche (absolute loss of 0.29 over the period 1970–2021).
• In 78 departments (over 96), the SWI trend between 1970–2021 is negative, and in 33 departments the trend is particularly severe. In these departments, the absolute SWI loss is higher than 0.1 (all seasons averaged).

Complementary notes — contact

This article and the code has been elaborated with the support of Ageoce. Ageoce builds digital solutions and services focusing on geodata intelligence and geosciences.

Please feel free to contact us anytime.

Appendix 1 : mathematical definition of the SWI

The calculation of the SWI reads as follow:

SWI = (θfc ​− θwi) / (​θ − θwi)​​

with:

• θ: water content in the soil,
• θwi: water content in the soil at the wilting point. The wilting point represents the point below what water cannot be extracted by plant roots,
• θfc water content in the soil at the field capacity. The field capacity represents the point after which water cannot be stored by soil any more. After that point, gravitational forces become too high and water starts to infiltrate the lower levels.

Additionally, for a month m, the uniformized SWI is given as follow:

USWI = (SWIm ​+ SWIm−1​ + SWIm−2) / 3​​

Appendix 2: Table of SWI evolution by department

Bibliography

[1] SDES (2021) Nouveau zonage d’exposition au retrait-gonflement des argiles : plus de 10,4 millions de maisons individuelles potentiellement très exposées. Article web du ministère en charge de l’écologie en France.

[2] France Assureurs (2022) Le risque sécheresse et son impact sur les habitations. Web article.

[3] Ministère de l’intérieur (2019) Circulaire INTE1911312C relative à la procédure de reconnaissance de l’état de catastrophe naturelle — Révision des critères permettant de caractériser l’intensité des épisodes de sécheresse-réhydratation des sols à l’origine de mouvements de terrains différentiels.

Groundwater recharge represents the amount of water coming from precipitation reaching the groundwater table. Its determination helps to better understand the available/renewable groundwater in watersheds and the shape of groundwater flow systems.

One of the simplest methods to estimate groundwater recharge is the Thornthwaite-Mather procedure (Steenhuis and Van Der Molen, 1986). This procedure was published by Thornthwaite and Mather (1955, 1957). The idea of this procedure is to calculate the water balance in the root zone of the soil where water can be (1) evaporated into the atmosphere under the effect of heat, (2) transpired by vegetation, (3) stored by the soil, and eventually (4) infiltrated when stored water exceeds the field capacity.

This procedures relies on several parameters and variables described as follows:

• information about soil texture (e.g. sand and clay content) to describe the hydraulic properties of the soil and its capacity to store/infiltrate,
• meteorological records: precipitation and potential evapotranspiration.

Of course groundwater recharge can be influenced by many other factors such as the slope of the terrain, the snow cover, the variability of the crop/land cover and the irrigation. In the following these aspects are not taken into account.

In the first part of the tutorial, the Earth Engine python API will be initialized, some useful libraries will be imported, and the location/period of interest will be defined. If you are not familiar with the Earth Engine python API, having a first look at the Intro to the Earth Engine python API is recommanded.

In the second part, OpenLandMap datasets related to soil properties will be explored. The wilting point and field capacity of the soil will be calculated by applying some mathematical expressions to multiple images.

In the third part, evapotranspiration and precipitation datasets will be imported. A function will be defined to resample the time resolution of an ee.ImageCollection and to homogenize time index of both datasets. Both datasets will then be combined into one.

In the fourth and final part, the Thornthwaite-Mather(TM) procedure will be implemented by iterating over the meteorological ee.ImageCollection. Finally, a comparison between groundwater recharge in two places will be described and the resulting mean annual groundwater recharge will be displayed over France.

Run me first

Earth Engine API

First of all, run the following cell to initialize the API. The output will contain instructions on how to grant this notebook access to Earth Engine using your account.

import ee

# Trigger the authentication flow.
ee.Authenticate()

# Initialize the library.
ee.Initialize()

Other libraries

Import other libraries/modules used in this notebook.

• pandas: data analysis (including the DataFrame data structure)
• matplotlib: data visualization library
• numpy: array-processing package
• folium: interactive web map
• pprint: a pretty printer
• branca.colormap: utility module for dealing with colormaps.

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import folium
import pprint
import branca.colormap as cm

Some input parameters

We additionally define some parameters used to evaluate the results of our implementation. Particularly:

• the period of interest to get meteorological records,
• a location of interest based on longitude and latitude coordinates. In the following, the point of interest is located in a productive agricultural region which is about 30 kilometers outside of the city of Lyon (France). This point is used to evaluate and illustrate the progress of the described procedure.

# Initial date of interest (inclusive).
i_date = "2015-01-01"

# Final date of interest (exclusive).
f_date = "2020-01-01"

# Define the location of interest with a point.
lon = 5.145041
lat = 45.772439
poi = ee.Geometry.Point(lon, lat)

# A nominal scale in meters of the projection to work in [in meters].
scale = 1000

From soil texture to hydraulic properties

Definitions

Two hydraulic properties of soil are commonly used in the TM procedure:

• the wilting point represents the point below what water cannot be extracted by plant roots,
• the field capacity represents the point after which water cannot be stored by soil any more. After that point, gravitational forces become too high and water starts to infiltrate the lower levels.

Some equations given by Saxton & Rawls (2006) are used to link both parameters to the texture of the soil. The calculation of water content at wilting point θWP can be done as follows:

θWP = θ1500t + (0.14 θ1500t − 0.002) with:

θ1500t = −0.024 S + 0.487 C + 0.006 OM + 0.005(S × OM)− 0.013(C × OM) + 0.068(S × C) + 0.031

where:

• S: represents the sand content of the soil (mass percentage),
• C: represents the clay content of the soil (mass percentage),
• OM: represents the organic matter content of the soil (mass percentage).

Similarly, the calculation of the water content at field capacity θFC can be done as follows:

θFC = θ33t + (1.283 θ33t² − 0.374 θ33t −0.15) with:

θ33t = -0.251 S + 0.195 C + 0.011 OM + 0.006 (S × OM) — 0.027 (C × OM) + 0.452 (S × C) + 0.299

Determination of soil texture and properties

In the following, OpenLandMap datasets are used to describe clay, sand and organic carbon content of soil. A global dataset of soil water content at the field capacity with a resolution of 250 m has been made available by Hengl & Gupta (2019). However, up to now, there is no dataset dedicated to the water content of soil at the wilting point. Consequently, in the following, both parameters will be determined considering the previous equations and using the global datasets giving the sand, clay and organic matter contents of the soil. According to the description, these datasets are based on machine learning predictions from global compilation of soil profiles and samples. Processing steps are described in detail here. The information (clay, sand content, etc.) is given at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution.

These standard depths and associated bands are defined into a list as follows:

# Soil depths [in cm] where we have data.
olm_depths = [0, 10, 30, 60, 100, 200]

# Names of bands associated with reference depths.
olm_bands = ["b" + str(sd) for sd in olm_depths]

We now define a function to get the ee.Image associated to the parameter we are interested in (e.g. sand, clay, organic carbon content, etc.).

def get_soil_prop(param):
    """
    This function returns soil properties image
    param (str): must be one of:
        "sand"     - Sand fraction
        "clay"     - Clay fraction
        "orgc"     - Organic Carbon fraction
    """
    if param == "sand":  # Sand fraction [%w]
        snippet = "OpenLandMap/SOL/SOL_SAND-WFRACTION_USDA-3A1A1A_M/v02"
        # Define the scale factor in accordance with the dataset description.
        scale_factor = 1 * 0.01

    elif param == "clay":  # Clay fraction [%w]
        snippet = "OpenLandMap/SOL/SOL_CLAY-WFRACTION_USDA-3A1A1A_M/v02"
        # Define the scale factor in accordance with the dataset description.
        scale_factor = 1 * 0.01

    elif param == "orgc":  # Organic Carbon fraction [g/kg]
        snippet = "OpenLandMap/SOL/SOL_ORGANIC-CARBON_USDA-6A1C_M/v02"
        # Define the scale factor in accordance with the dataset description.
        scale_factor = 5 * 0.001  # to get kg/kg
    else:
        return print("error")

    # Apply the scale factor to the ee.Image.
    dataset = ee.Image(snippet).multiply(scale_factor)

    return dataset

We apply this function to import soil properties:

# Image associated with the sand content.
sand = get_soil_prop("sand")

# Image associated with the clay content.
clay = get_soil_prop("clay")

# Image associated with the organic carbon content.
orgc = get_soil_prop("orgc")

To illustrate the result, we define a new method for handing Earth Engine tiles and using it to display the clay content of the soil at a given reference depth, to a Leaflet map.

def add_ee_layer(self, ee_image_object, vis_params, name):
    """Adds a method for displaying Earth Engine image tiles to folium map."""
    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)
    folium.raster_layers.TileLayer(
        tiles=map_id_dict["tile_fetcher"].url_format,
        attr="Map Data &copy; <a href='https://earthengine.google.com/'>Google Earth Engine</a>",
        name=name,
        overlay=True,
        control=True,
    ).add_to(self)


# Add Earth Engine drawing method to folium.
folium.Map.add_ee_layer = add_ee_layer

my_map = folium.Map(location=[lat, lon], zoom_start=3)

# Set visualization parameters.
vis_params = {
    "bands": ["b0"],
    "min": 0.01,
    "max": 1,
    "opacity": 1,
    "palette": ["white", "#464646"],
}

# Add the sand content data to the map object.
my_map.add_ee_layer(sand, vis_params, "Sand Content")

# Add a marker at the location of interest.
folium.Marker([lat, lon], popup="point of interest").add_to(my_map)

# Add a layer control panel to the map.
my_map.add_child(folium.LayerControl())

# Display the map.
display(my_map)

Now, a function is defined to get soil properties at a given location. The following function returns a dictionary indicating the value of the parameter of interest for each standard depth (in centimeter). This function uses the ee.Image.sample method to evaluate the ee.Image properties on the region of interest. The result is then transferred client-side using the ee.Image.getInfo method.

In the example below, we are asking for the sand content.

def local_profile(dataset, poi, buffer):
    # Get properties at the location of interest and transfer to client-side.
    prop = dataset.sample(poi, buffer).select(olm_bands).getInfo()

    # Selection of the features/properties of interest.
    profile = prop["features"][0]["properties"]

    # Re-shaping of the dict.
    profile = {key: round(val, 3) for key, val in profile.items()}

    return profile


# Apply the function to get the sand profile.
profile_sand = local_profile(sand, poi, scale)

# Print the result.
print("Sand content profile at the location of interest:\n", profile_sand)

We now apply the function to plot the profile of the soil regarding sand and clay and organic carbon content at the location of interest:

# Clay and organic content profiles.
profile_clay = local_profile(clay, poi, scale)
profile_orgc = local_profile(orgc, poi, scale)

# Data visualization in the form of a bar plot.
fig, ax = plt.subplots(figsize=(15, 6))
ax.axes.get_yaxis().set_visible(False)

# Definition of label locations.
x = np.arange(len(olm_bands))

# Definition of the bar width.
width = 0.25

# Bar plot representing the sand content profile.
rect1 = ax.bar(
    x - width,
    [round(100 * profile_sand[b], 2) for b in olm_bands],
    width,
    label="Sand",
    color="#ecebbd",
)

# Bar plot representing the clay content profile.
rect2 = ax.bar(
    x,
    [round(100 * profile_clay[b], 2) for b in olm_bands],
    width,
    label="Clay",
    color="#6f6c5d",
)

# Bar plot representing the organic carbon content profile.
rect3 = ax.bar(
    x + width,
    [round(100 * profile_orgc[b], 2) for b in olm_bands],
    width,
    label="Organic Carbon",
    color="black",
    alpha=0.75,
)

# Definition of a function to attach a label to each bar.
def autolabel_soil_prop(rects):
    """Attach a text label above each bar in *rects*, displaying its height."""
    for rect in rects:
        height = rect.get_height()
        ax.annotate(
            "{}".format(height) + "%",
            xy=(rect.get_x() + rect.get_width() / 2, height),
            xytext=(0, 3),  # 3 points vertical offset.
            textcoords="offset points",
            ha="center",
            va="bottom",
            fontsize=10,
        )


# Application of the function to each barplot.
autolabel_soil_prop(rect1)
autolabel_soil_prop(rect2)
autolabel_soil_prop(rect3)

# Title of the plot.
ax.set_title("Properties of the soil at different depths (mass content)", fontsize=14)

# Properties of x/y labels and ticks.
ax.set_xticks(x)
x_labels = [str(d) + " cm" for d in olm_depths]
ax.set_xticklabels(x_labels, rotation=45, fontsize=10)

ax.spines["left"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)

# Shrink current axis's height by 10% on the bottom.
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1, box.width, box.height * 0.9])

# Add a legend below current axis.
ax.legend(
    loc="upper center", bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=3
)

plt.show()

Expression to calculate hydraulic properties

Now that soil properties are described, the water content at the field capacity and at the wilting point can be calculated according to the equation defined at the beginning of this section. Please note that in the equation of Saxton & Rawls (2006), the wilting point and field capacity are calculated using the Organic Matter content (OM) and not the Organic Carbon content (OC). In the following, we convert OC into OM using the corrective factor known as the Van Bemmelen factor:

0M = 1.724 × OC

Several operators are available to perform basic mathematical operations on image bands: add(), subtract(), multiply() and divide(). Here, we multiply the organic content by the Van Bemmelen factor. It is done using the ee.Image.multiply method on the organic carbon content ee.Image.

# Conversion of organic carbon content into organic matter content.
orgm = orgc.multiply(1.724)

# Organic matter content profile.
profile_orgm = local_profile(orgm, poi, scale)
print("Organic Matter content profile at the location of interest:\n", profile_orgm)

When the mathematical operation to apply to the ee.Image becomes too complex, the ee.Image.expression is a good alternative. We use it in the following code block since the calculation of wilting point and field capacity relies on multiple parameters and images. This method takes two arguments:

• a string formalizing the arithmetic expression we want to evaluate,
• a dict associating images to each parameter of the arithmetic expression.

The mathematical expression is applied as follows to determine wilting point and field capacity:

# Initialization of two constant images for wilting point and field capacity.
wilting_point = ee.Image(0)
field_capacity = ee.Image(0)

# Calculation for each standard depth using a loop.
for key in olm_bands:
    # Getting sand, clay and organic matter at the appropriate depth.
    si = sand.select(key)
    ci = clay.select(key)
    oi = orgm.select(key)

    # Calculation of the wilting point.
    # The theta_1500t parameter is needed for the given depth.
    theta_1500ti = (
        ee.Image(0)
        .expression(
            "-0.024 * S + 0.487 * C + 0.006 * OM + 0.005 * (S * OM)\
        - 0.013 * (C * OM) + 0.068 * (S * C) + 0.031",
            {
                "S": si,
                "C": ci,
                "OM": oi,
            },
        )
        .rename("T1500ti")
    )

    # Final expression for the wilting point.
    wpi = theta_1500ti.expression(
        "T1500ti + ( 0.14 * T1500ti - 0.002)", {"T1500ti": theta_1500ti}
    ).rename("wpi")

    # Add as a new band of the global wilting point ee.Image.
    # Do not forget to cast the type with float().
    wilting_point = wilting_point.addBands(wpi.rename(key).float())

    # Same process for the calculation of the field capacity.
    # The parameter theta_33t is needed for the given depth.
    theta_33ti = (
        ee.Image(0)
        .expression(
            "-0.251 * S + 0.195 * C + 0.011 * OM +\
        0.006 * (S * OM) - 0.027 * (C * OM)+\
        0.452 * (S * C) + 0.299",
            {
                "S": si,
                "C": ci,
                "OM": oi,
            },
        )
        .rename("T33ti")
    )

    # Final expression for the field capacity of the soil.
    fci = theta_33ti.expression(
        "T33ti + (1.283 * T33ti * T33ti - 0.374 * T33ti - 0.015)",
        {"T33ti": theta_33ti.select("T33ti")},
    )

    # Add a new band of the global field capacity ee.Image.
    field_capacity = field_capacity.addBands(fci.rename(key).float())

Let’s see the result around our location of interest:

profile_wp = local_profile(wilting_point, poi, scale)
profile_fc = local_profile(field_capacity, poi, scale)

print("Wilting point profile:\n", profile_wp)
print("Field capacity profile:\n", profile_fc)

The result is displayed using barplots as follows:

fig, ax = plt.subplots(figsize=(15, 6))
ax.axes.get_yaxis().set_visible(False)

# Definition of the label locations.
x = np.arange(len(olm_bands))

# Width of the bar of the barplot.
width = 0.25

# Barplot associated with the water content at the wilting point.
rect1 = ax.bar(
    x - width / 2,
    [round(profile_wp[b] * 100, 2) for b in olm_bands],
    width,
    label="Water content at wilting point",
    color="red",
    alpha=0.5,
)

# Barplot associated with the water content at the field capacity.
rect2 = ax.bar(
    x + width / 2,
    [round(profile_fc[b] * 100, 2) for b in olm_bands],
    width,
    label="Water content at field capacity",
    color="blue",
    alpha=0.5,
)

# Add Labels on top of bars.
autolabel_soil_prop(rect1)
autolabel_soil_prop(rect2)

# Title of the plot.
ax.set_title("Hydraulic properties of the soil at different depths", fontsize=14)

# Properties of x/y labels and ticks.
ax.set_xticks(x)
x_labels = [str(d) + " cm" for d in olm_depths]
ax.set_xticklabels(x_labels, rotation=45, fontsize=10)

ax.spines["left"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)

# Shrink current axis's height by 10% on the bottom.
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1, box.width, box.height * 0.9])

# Put a legend below current axis.
ax.legend(
    loc="upper center", bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=2
)

plt.show()

Getting meteorological datasets

Datasets exploration

The meteorological data used in our implementation of the TM procedure relies on the following datasets:

• Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) gives precipitations on a daily basis (resolution of 5 km),
• MODIS Terra Net gives evapotranspiration on a 8-days basis (resolution of 500 m).

Both datasets are imported as follows, specifying the bands of interest using .select() and the period of interest using .filterDate().

# Import precipitation.
pr = (
    ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
    .select("precipitation")
    .filterDate(i_date, f_date)
)

# Import potential evaporation PET and its quality indicator ET_QC.
pet = (
    ee.ImageCollection("MODIS/006/MOD16A2")
    .select(["PET", "ET_QC"])
    .filterDate(i_date, f_date)
)

Now we can have a closer look around our location of interest. To evaluate the properties of an ee.ImageCollection, the ee.ImageCollection.getRegion method is used and combined with ee.ImageCollection.getInfo method for a client-side visualization.

# Evaluate local precipitation conditions.
local_pr = pr.getRegion(poi, scale).getInfo()
pprint.pprint(local_pr[:5])

We now establish a procedure to get meteorological data around a given location in form of a pandas.DataFrame:

def ee_array_to_df(arr, list_of_bands):
    """Transforms client-side ee.Image.getRegion array to pandas.DataFrame."""
    df = pd.DataFrame(arr)

    # Rearrange the header.
    headers = df.iloc[0]
    df = pd.DataFrame(df.values[1:], columns=headers)

    # Convert the data to numeric values.
    for band in list_of_bands:
        df[band] = pd.to_numeric(df[band], errors="coerce")

    # Convert the time field into a datetime.
    df["datetime"] = pd.to_datetime(df["time"], unit="ms")

    # Keep the columns of interest.
    df = df[["time", "datetime", *list_of_bands]]

    # The datetime column is defined as index.
    df = df.set_index("datetime")

    return df

We apply the function and see the head of the resulting pandas.DataFrame:

pr_df = ee_array_to_df(local_pr, ["precipitation"])
pr_df.head(10)

We do the same for potential evaporation:

# Evaluate local potential evapotranspiration.
local_pet = pet.getRegion(poi, scale).getInfo()

# Transform the result into a pandas dataframe.
pet_df = ee_array_to_df(local_pet, ["PET", "ET_QC"])
pet_df.head(5)

Looking at both pandas.DataFrame shows the following points:

• the time resolution between both datasets is not the same,
• for some reasons, potential evapotranspiration cannot be calculated at some dates. It corresponds to lines where the quality indicator ET_QC is higher than 1.

Both issues must be handled before implementing the iterative process: we want to work on a similar timeline with potential evapotranspiration and precipitation, and we want to avoid missing values.

Resampling the time resolution of an ee.ImageCollection

To address these issues (homogeneous time index and missing values), we make a sum resampling of both datasets by month. When PET cannot be calculated, the monthly averaged value is considered. The key steps and functions used to resample are described below:

• A new date index is defined as a sequence using the ee.List.sequence method.
• A function representing the resampling operation is defined. This function consists of grouping images of the desired time interval and calculating the sum. The sum is calculated by taking the mean between available images and multiplying it by the duration of the interval.
• The user-supplied function is then mapped over the new time index using .map().

Finally, the resampling procedure reads as follows:

def sum_resampler(coll, freq, unit, scale_factor, band_name):
    """
    This function aims to resample the time scale of an ee.ImageCollection.
    The function returns an ee.ImageCollection with the averaged sum of the
    band on the selected frequency.

    coll: (ee.ImageCollection) only one band can be handled
    freq: (int) corresponds to the resampling frequence
    unit: (str) corresponds to the resampling time unit.
                must be 'day', 'month' or 'year'
    scale_factor (float): scaling factor used to get our value in the good unit
    band_name (str) name of the output band
    """

    # Define initial and final dates of the collection.
    firstdate = ee.Date(
        coll.sort("system:time_start", True).first().get("system:time_start")
    )

    lastdate = ee.Date(
        coll.sort("system:time_start", False).first().get("system:time_start")
    )

    # Calculate the time difference between both dates.
    # https://developers.google.com/earth-engine/apidocs/ee-date-difference
    diff_dates = lastdate.difference(firstdate, unit)

    # Define a new time index (for output).
    new_index = ee.List.sequence(0, ee.Number(diff_dates), freq)

    # Define the function that will be applied to our new time index.
    def apply_resampling(date_index):
        # Define the starting date to take into account.
        startdate = firstdate.advance(ee.Number(date_index), unit)

        # Define the ending date to take into account according
        # to the desired frequency.
        enddate = firstdate.advance(ee.Number(date_index).add(freq), unit)

        # Calculate the number of days between starting and ending days.
        diff_days = enddate.difference(startdate, "day")

        # Calculate the composite image.
        image = (
            coll.filterDate(startdate, enddate)
            .mean()
            .multiply(diff_days)
            .multiply(scale_factor)
            .rename(band_name)
        )

        # Return the final image with the appropriate time index.
        return image.set("system:time_start", startdate.millis())

    # Map the function to the new time index.
    res = new_index.map(apply_resampling)

    # Transform the result into an ee.ImageCollection.
    res = ee.ImageCollection(res)

    return res

The precipitation dataset is now resampled by month as follows:

• the collection to resample is defined as pr,
• we want a collection on a monthly basis then freq = 1 and unit = "month",
• there is no correction factor to apply according to the dataset description then scale_factor = 1,
• "pr" is the name of the output band.

# Apply the resampling function to the precipitation dataset.
pr_m = sum_resampler(pr, 1, "month", 1, "pr")

# Evaluate the result at the location of interest.
pprint.pprint(pr_m.getRegion(poi, scale).getInfo()[:5])

For evapotranspiration, we have to be careful with the unit. The dataset gives us an 8-day sum and a scale factor of 10 is applied. Then, to get a homogeneous unit, we need to rescale by dividing by 8 and 10 = 0.0125.

# Apply the resampling function to the PET dataset.
pet_m = sum_resampler(pet.select("PET"), 1, "month", 0.0125, "pet")

# Evaluate the result at the location of interest.
pprint.pprint(pet_m.getRegion(poi, scale).getInfo()[:5])

We now combine both ee.ImageCollection objects (pet_m and pr_m) using the ee.ImageCollection.combine method. Note that corresponding images in both ee.ImageCollection objects need to have the same time index before combining.

# Combine precipitation and evapotranspiration.
meteo = pr_m.combine(pet_m)

# Import meteorological data as an array at the location of interest.
meteo_arr = meteo.getRegion(poi, scale).getInfo()

# Print the result.
pprint.pprint(meteo_arr[:5])

We evaluate the result on our location of interest:

# Transform the array into a pandas dataframe and sort the index.
meteo_df = ee_array_to_df(meteo_arr, ["pr", "pet"]).sort_index()

# Data visualization
fig, ax = plt.subplots(figsize=(15, 6))

# Barplot associated with precipitations.
meteo_df["pr"].plot(kind="bar", ax=ax, label="precipitation")

# Barplot associated with potential evapotranspiration.
meteo_df["pet"].plot(
    kind="bar", ax=ax, label="potential evapotranspiration", color="orange", alpha=0.5
)

# Add a legend.
ax.legend()

# Add some x/y-labels properties.
ax.set_ylabel("Intensity [mm]")
ax.set_xlabel(None)

# Define the date format and shape of x-labels.
x_labels = meteo_df.index.strftime("%m-%Y")
ax.set_xticklabels(x_labels, rotation=90, fontsize=10)

plt.show()

Implementation of the TM procedure

Description

Some additional definitions are needed to formalize the Thornthwaite-Mather procedure. The following definitions are given in accordance with Allen et al. (1998) (the document can be downloaded here):

TAW = 1000 × (θFC − θWP) × Zr where:

• TAW: the total available soil water in the root zone [mm],
• θFC: the water content at the field capacity [m³ m³],
• θWP: the water content at wilting point [m³ m³],
• Zr: the rooting depth [mm],

Typical values of θFC and θWP for different soil types are given in Table 19 of Allen et al. (1998).

The readily available water (RAW) is given by RAW = p × TAW, where p is the average fraction of TAW that can be depleted from the root zone before moisture stress (ranging between 0 to 1). This quantity is also noted STFC which is the available water stored at field capacity in the root zone.

Ranges of maximum effective rooting depth Zr, and soil water depletion fraction for no stress p, for common crops are given in the Table 22 of Allen et al. (1998). In addition, a global effective plant rooting depth dataset is provided by Yang et al. (2016) with a resolution of 0.5° by 0.5° (see the paper here and the dataset here).

According to this global dataset, the effective rooting depth around our region of interest (France) can reasonably assumed to Zr = 0.5. Additionally, the parameter p is also assumed constant and equal to and p = 0.5 which is in line with common values described in Table 22 of Allen et al. (1998).

zr = ee.Image(0.5)
p = ee.Image(0.5)

In the following, we also consider an averaged value between reference depths of the water content at wilting point and field capacity:

def olm_prop_mean(olm_image, band_output_name):
    """
    This function calculates an averaged value of
    soil properties between reference depths.
    """
    mean_image = olm_image.expression(
        "(b0 + b10 + b30 + b60 + b100 + b200) / 6",
        {
            "b0": olm_image.select("b0"),
            "b10": olm_image.select("b10"),
            "b30": olm_image.select("b30"),
            "b60": olm_image.select("b60"),
            "b100": olm_image.select("b100"),
            "b200": olm_image.select("b200"),
        },
    ).rename(band_output_name)

    return mean_image


# Apply the function to field capacity and wilting point.
fcm = olm_prop_mean(field_capacity, "fc_mean")
wpm = olm_prop_mean(wilting_point, "wp_mean")

# Calculate the theoretical available water.
taw = (
    (fcm.select("fc_mean").subtract(wpm.select("wp_mean"))).multiply(1000).multiply(zr)
)

# Calculate the stored water at the field capacity.
stfc = taw.multiply(p)

The Thornthwaite-Mather procedure used to estimate groundwater recharge is explicitly described by Steenhuis and Van der Molen (1985). This procedure uses monthly sums of potential evaporation, cumulative precipitation, and the moisture status of the soil which is calculated iteratively. The moisture status of the soils depends on the accumulated potential water loss (APWL). This parameter is calculated depending on whether the potential evaporation is greater than or less than the cumulative precipitation. The procedure reads as follow:

Case 1: potential evapotranspiration is higher than precipitation.

In that case, PET>P and APWLm is incremented as follows: APWLm = APWLm-1+ (PETm — Pm) where:

• APWLm (respectively APWLm-1) represents the accumulated potential water loss for the month m (respectively at the previous month m-1)
• PETm the cumulative potential evapotranspiration at month m,
• Pm the cumulative precipitation at month m,

and the relationship between APWL and the amount of water stored in the root zone for the month m is expressed as: STm = STFC × [exp(-APWLm/STFC)] where STm is the available water stored in the root zone for the month m.

Case 2: potential evapotranspiration is lower than precipitation.

In that case, PET<P and STm is incremented as follows:

STm = STm-1 + (Pm-PETm)

Case 2.1: the storage STm is higher than the water stored at the field capacity.

If STm > STFC the recharge is calculated as:

Rm = STm - STFC + Pm - PETm

In addition, the water stored at the end of the month m becomes equal to STFC and APWLm is set equal to zero.

Case 2.2: the storage STm is less than or equal to the water stored at the field capacity.

If STm < STFC, APWLm is implemented as follows:

APWLm = STFC × ln(STm/STFC), and no percolation occurs.

Initialization

The initial time of the calculation is defined according to the first date of the meteorological collection:

# Define the initial time (time0) according to the start of the collection.
time0 = meteo.first().get("system:time_start")

Then, we initialize the calculation with an ee.Image where all bands associated to the hydric state of the soil are set equal to ee.Image(0), except for the initial storage which is considered to be equal to the water content at field capacity, meaning that ST0 = STFC.

# Initialize all bands describing the hydric state of the soil.
# Do not forget to cast the type of the data with a .float().
# Initial recharge.
initial_rech = ee.Image(0).set("system:time_start", time0).select([0], ["rech"]).float()

# Initialization of APWL.
initial_apwl = ee.Image(0).set("system:time_start", time0).select([0], ["apwl"]).float()

# Initialization of ST.
initial_st = stfc.set("system:time_start", time0).select([0], ["st"]).float()

# Initialization of precipitation.
initial_pr = ee.Image(0).set("system:time_start", time0).select([0], ["pr"]).float()

# Initialization of potential evapotranspiration.
initial_pet = ee.Image(0).set("system:time_start", time0).select([0], ["pet"]).float()

We combine all these bands into one ee.Image adding new bands to the first using the ee.Image.addBands method:

initial_image = initial_rech.addBands(
    ee.Image([initial_apwl, initial_st, initial_pr, initial_pet])
)

We also initialize a list in which new images will be added after each iteration. We create this server-side list using the ee.List method.

image_list = ee.List([initial_image])

Iteration over an ee.ImageCollection

The procedure is implemented by means of the ee.ImageCollection.iterate method, which applies a user-supplied function to each element of a collection. For each time step, groundwater recharge is calculated using the recharge_calculator considering the previous hydric state of the soil and current meteorological conditions.

Of course, considering the TM description, several cases must be distinguished to calculate groundwater recharge. The distinction is made by the definition of binary layers with different logical operations. It allows specific calculations to be applied in areas where a given condition is true using the ee.Image.where method.

The function we apply to each element of the meteorological dataset to calculate groundwater recharge is defined as follows.

def recharge_calculator(image, image_list):
    """
    Contains operations made at each iteration.
    """
    # Determine the date of the current ee.Image of the collection.
    localdate = image.date().millis()

    # Import previous image stored in the list.
    prev_im = ee.Image(ee.List(image_list).get(-1))

    # Import previous APWL and ST.
    prev_apwl = prev_im.select("apwl")
    prev_st = prev_im.select("st")

    # Import current precipitation and evapotranspiration.
    pr_im = image.select("pr")
    pet_im = image.select("pet")

    # Initialize the new bands associated with recharge, apwl and st.
    # DO NOT FORGET TO CAST THE TYPE WITH .float().
    new_rech = (
        ee.Image(0)
        .set("system:time_start", localdate)
        .select([0], ["rech"])
        .float()
    )

    new_apwl = (
        ee.Image(0)
        .set("system:time_start", localdate)
        .select([0], ["apwl"])
        .float()
    )

    new_st = (
        prev_st.set("system:time_start", localdate).select([0], ["st"]).float()
    )

    # Calculate bands depending on the situation using binary layers with
    # logical operations.

    # CASE 1.
    # Define zone1: the area where PET > P.
    zone1 = pet_im.gt(pr_im)

    # Calculation of APWL in zone 1.
    zone1_apwl = prev_apwl.add(pet_im.subtract(pr_im)).rename("apwl")
    # Implementation of zone 1 values for APWL.
    new_apwl = new_apwl.where(zone1, zone1_apwl)

    # Calculate ST in zone 1.
    zone1_st = prev_st.multiply(
        ee.Image.exp(zone1_apwl.divide(stfc).multiply(-1))
    ).rename("st")
    # Implement ST in zone 1.
    new_st = new_st.where(zone1, zone1_st)

    # CASE 2.
    # Define zone2: the area where PET <= P.
    zone2 = pet_im.lte(pr_im)

    # Calculate ST in zone 2.
    zone2_st = prev_st.add(pr_im).subtract(pet_im).rename("st")
    # Implement ST in zone 2.
    new_st = new_st.where(zone2, zone2_st)

    # CASE 2.1.
    # Define zone21: the area where PET <= P and ST >= STfc.
    zone21 = zone2.And(zone2_st.gte(stfc))

    # Calculate recharge in zone 21.
    zone21_re = zone2_st.subtract(stfc).rename("rech")
    # Implement recharge in zone 21.
    new_rech = new_rech.where(zone21, zone21_re)
    # Implement ST in zone 21.
    new_st = new_st.where(zone21, stfc)

    # CASE 2.2.
    # Define zone 22: the area where PET <= P and ST < STfc.
    zone22 = zone2.And(zone2_st.lt(stfc))

    # Calculate APWL in zone 22.
    zone22_apwl = (
        stfc.multiply(-1).multiply(ee.Image.log(zone2_st.divide(stfc))).rename("apwl")
    )
    # Implement APWL in zone 22.
    new_apwl = new_apwl.where(zone22, zone22_apwl)

    # Create a mask around area where recharge can effectively be calculated.
    # Where we have have PET, P, FCm, WPm (except urban areas, etc.).
    mask = pet_im.gte(0).And(pr_im.gte(0)).And(fcm.gte(0)).And(wpm.gte(0))

    # Apply the mask.
    new_rech = new_rech.updateMask(mask)

    # Add all Bands to our ee.Image.
    new_image = new_rech.addBands(ee.Image([new_apwl, new_st, pr_im, pet_im]))

    # Add the new ee.Image to the ee.List.
    return ee.List(image_list).add(new_image)

The TM procedure can now be applied to the meteorological ee.ImageCollection:

# Iterate the user-supplied function to the meteo collection.
rech_list = meteo.iterate(recharge_calculator, image_list)

# Remove the initial image from our list.
rech_list = ee.List(rech_list).remove(initial_image)

# Transform the list into an ee.ImageCollection.
rech_coll = ee.ImageCollection(rech_list)

Let’s have a look at the result around the location of interest:

arr = rech_coll.getRegion(poi, scale).getInfo()
rdf = ee_array_to_df(arr, ["pr", "pet", "apwl", "st", "rech"]).sort_index()
rdf.head(12)

The result can be displayed in the form of a barplot as follows:

# Data visualization in the form of barplots.
fig, ax = plt.subplots(figsize=(15, 6))

# Barplot associated with precipitation.
rdf["pr"].plot(kind="bar", ax=ax, label="precipitation", alpha=0.5)

# Barplot associated with potential evapotranspiration.
rdf["pet"].plot(
    kind="bar", ax=ax, label="potential evapotranspiration", color="orange", alpha=0.2
)

# Barplot associated with groundwater recharge
rdf["rech"].plot(kind="bar", ax=ax, label="recharge", color="green", alpha=1)

# Add a legend.
ax.legend()

# Define x/y-labels properties.
ax.set_ylabel("Intensity [mm]")
ax.set_xlabel(None)

# Define the date format and shape of x-labels.
x_labels = rdf.index.strftime("%m-%Y")
ax.set_xticklabels(x_labels, rotation=90, fontsize=10)

plt.show()

The result shows the distribution of precipitation, potential evapotranspiration, and groundwater recharge along the year. It shows that in our area of interest, groundwater recharge generally occurs from October to March. Even though a significant amount of precipitation occurs from April to September, evapotranspiration largely dominates because of high temperatures and sun exposure during these months. The result is that no percolation into aquifers occurs during this period.

Now the annual average recharge over the period of interest can be calculated. To do that, we resample the DataFrame we've just created:

# Resample the pandas dataframe on a yearly basis making the sum by year.
rdfy = rdf.resample("Y").sum()

# Calculate the mean value.
mean_recharge = rdfy["rech"].mean()

# Print the result.
print(
    "The mean annual recharge at our point of interest is", int(mean_recharge), "mm/an"
)

The mean annual recharge at our point of interest is 154 mm/an

Groundwater recharge comparison between multiple places

We now may want to get local information about groundwater recharge and/or map this variable on an area of interest.

Let’s define a function to get the local information based on the ee.ImageCollection we've just built:

def get_local_recharge(i_date, f_date, lon, lat, scale):
    """
    Returns a pandas df describing the cumulative groundwater
    recharge by month
    """
    # Define the location of interest with a point.
    poi = ee.Geometry.Point(lon, lat)

    # Evaluate the recharge around the location of interest.
    rarr = rech_coll.filterDate(i_date, f_date).getRegion(poi, scale).getInfo()

    # Transform the result into a pandas dataframe.
    rdf = ee_array_to_df(rarr, ["pr", "pet", "apwl", "st", "rech"]).sort_index()
    return rdf

We now use this function on a second point of interest located near the city of Montpellier (France). This city is located in the south of France, and precipitation and groundwater recharge are expected to be much lower than in the previous case.

# Define the second location of interest by longitude/latitude.
lon2 = 4.137152
lat2 = 43.626945

# Calculate the local recharge condition at this location.
rdf2 = get_local_recharge(i_date, f_date, lon2, lat2, scale)

# Resample the resulting pandas dataframe on a yearly basis (sum by year).
rdf2y = rdf2.resample("Y").sum()
rdf2y.head()

# Data Visualization
fig, ax = plt.subplots(figsize=(15, 6))
ax.axes.get_yaxis().set_visible(False)

# Define the x-label locations.
x = np.arange(len(rdfy))

# Define the bar width.
width = 0.25

# Bar plot associated to groundwater recharge at the 1st location of interest.
rect1 = ax.bar(
    x - width / 2, rdfy.rech, width, label="Lyon (France)", color="blue", alpha=0.5
)

# Bar plot associated to groundwater recharge at the 2nd location of interest.
rect2 = ax.bar(
    x + width / 2,
    rdf2y.rech,
    width,
    label="Montpellier (France)",
    color="red",
    alpha=0.5,
)

# Define a function to attach a label to each bar.
def autolabel_recharge(rects):
    """Attach a text label above each bar in *rects*, displaying its height."""
    for rect in rects:
        height = rect.get_height()
        ax.annotate(
            "{}".format(int(height)) + " mm",
            xy=(rect.get_x() + rect.get_width() / 2, height),
            xytext=(0, 3),  # 3 points vertical offset
            textcoords="offset points",
            ha="center",
            va="bottom",
            fontsize=8,
        )


autolabel_recharge(rect1)
autolabel_recharge(rect2)

# Calculate the averaged annual recharge at both locations of interest.
place1mean = int(rdfy["rech"].mean())
place2mean = int(rdf2y["rech"].mean())

# Add an horizontal line associated with averaged annual values (location 1).
ax.hlines(
    place1mean,
    xmin=min(x) - width,
    xmax=max(x) + width,
    color="blue",
    lw=0.5,
    label="average " + str(place1mean) + " mm/y",
    alpha=0.5,
)

# Add an horizontal line associated with averaged annual values (location 2).
ax.hlines(
    place2mean,
    xmin=min(x) - width,
    xmax=max(x) + width,
    color="red",
    lw=0.5,
    label="average " + str(place2mean) + " mm/y",
    alpha=0.5,
)

# Add a title.
ax.set_title("Groundwater recharge comparison between two places", fontsize=12)

# Define some x/y-axis properties.
ax.set_xticks(x)
x_labels = rdfy.index.year.tolist()
ax.set_xticklabels(x_labels, rotation=45, fontsize=10)

ax.spines["left"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)

# Shrink current axis's height by 10% on the bottom.
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1, box.width, box.height * 0.9])

# Add a legend below current axis.
ax.legend(
    loc="upper center", bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=2
)

plt.show()

The result shows that the annual recharge in Lyon is almost twice as high as in the area of Montpellier. The result also shows a great variability of the annual recharge ranging from 98 mm/y to 258 mm/y in Lyon and from 16 mm/y to 147 mm/y in Montpellier.

Groundwater recharge map of France

To get a map of groundwater recharge around our region of interest, let’s create a mean composite ee.Image based on our resulting ee.ImageCollection.

# Calculate the averaged annual recharge.
annual_rech = rech_coll.select("rech").mean().multiply(12)

# Calculate the average annual precipitation.
annual_pr = rech_coll.select("pr").mean().multiply(12)

# Get a feature collection of administrative boundaries.
countries = ee.FeatureCollection("FAO/GAUL/2015/level0").select("ADM0_NAME")

# Filter the feature collection to subset France.
france = countries.filter(ee.Filter.eq("ADM0_NAME", "France"))

# Clip the composite ee.Images around the region of interest.
rech_france = annual_rech.clip(france)
pr_france = annual_pr.clip(france)

And finally the map can be drawn.

# Create a folium map.
my_map = folium.Map(location=[lat, lon], zoom_start=6, zoom_control=False)

# Set visualization parameters for recharge.
rech_vis_params = {
    "bands": "rech",
    "min": 0,
    "max": 500,
    "opacity": 1,
    "palette": ["red", "orange", "yellow", "green", "blue", "purple"],
}

# Set visualization parameters for precipitation.
pr_vis_params = {
    "bands": "pr",
    "min": 500,
    "max": 1500,
    "opacity": 1,
    "palette": ["white", "blue"],
}

# Define a recharge colormap.
rech_colormap = cm.LinearColormap(
    colors=rech_vis_params["palette"],
    vmin=rech_vis_params["min"],
    vmax=rech_vis_params["max"],
)

# Define a precipitation colormap.
pr_colormap = cm.LinearColormap(
    colors=pr_vis_params["palette"],
    vmin=pr_vis_params["min"],
    vmax=pr_vis_params["max"],
)

# Caption of the recharge colormap.
rech_colormap.caption = "Average annual recharge rate (mm/year)"

# Caption of the precipitation colormap.
pr_colormap.caption = "Average annual precipitation rate (mm/year)"

# Add the precipitation composite to the map object.
my_map.add_ee_layer(pr_france, pr_vis_params, "Precipitation")

# Add the recharge composite to the map object.
my_map.add_ee_layer(rech_france, rech_vis_params, "Recharge")

# Add a marker at both locations of interest.
folium.Marker([lat, lon], popup="Area of Lyon").add_to(my_map)
folium.Marker([lat2, lon2], popup="Area of Montpellier").add_to(my_map)

# Add the colormaps to the map.
my_map.add_child(rech_colormap)
my_map.add_child(pr_colormap)

# Add a layer control panel to the map.
my_map.add_child(folium.LayerControl())

# Display the map.
display(my_map)

The resulting map represents the averaged annual precipitation and groundwater recharge rate in France over the period 2015–2019. Groundwater recharge is represented between 0 (red) to 500 mm/year (purple). Precipitation is represented between 500 (white) to 1500 mm/year (blue).

Please note the following points:

• The precipitation dataset used in this tutorial is not available over +/- 50 deg latitudes. That is why (1) the northern part of France is not covered by the precipitation/recharge, and (2) some anomalies/discontinuities can be observed near the boundary of the dataset.
• The procedure described in this tutorial does not consider the terrain slope which has a strong influence on runoff, especially in areas with significant relief. In these areas, groundwater recharge is overestimated because runoff is not considered.
• Snow cover has a more complex contribution to groundwater recharge than precipitation. This process is not taken into account in this tutorial. Hence, groundwater recharge is overestimated in areas where snow cover period is significant.
• The recharge calculation can be highly influenced by missing precipitation or potential evapotranspiration values. This issue is partially addressed by the dataset resampling procedure, but missing value can still exist and this point must be checked before considering that the calculated recharge is representative and reliable.

Bibliography

Allen RG, Pereira LS, Raes D, Smith M (1998). Crop evapotranspiration: guidelines for computing crop water requirements. Irrigation and Drainage Paper 56, FAO, Rome.

Saxton, K. E., & Rawls, W. J. (2006). Soil water characteristic estimates by texture and organic matter for hydrologic solutions. Soil science society of America Journal, 70(5), 1569–1578.

Steenhuis, T. S., & Van der Molen, W. H. (1986). The Thornthwaite-Mather procedure as a simple engineering method to predict recharge. Journal of Hydrology, 84(3–4), 221–229.

Thornthwaite, C. W., & Mather, J. R. (1957). Instructions and tables for computing potential evapotranspiration and the water balance. Publ. Climatol., 10(3).

Yang, Y., Donohue, R. J., & McVicar, T. R. (2016). Global estimation of effective plant rooting depth: Implications for hydrological modeling. Water Resources Research, 52(10), 8260–8276.

Acknowledgements

Thanks to Susanne Benz for precious advice. Thanks to Justin Braaten and Gino Miceli for reviewing and helping write and publish this tutorial on the Earth Engine Community platform.

Complementary notes — contact

This article and the code has been elaborated with the support of Ageoce. Ageoce builds digital solutions and services focusing on geodata intelligence and geosciences.

Please feel free to contact us anytime.

Introduction

A Digital Elevation Model (DEM) is a description of the ground elevation. DEMs are used in many fields such as hydrology, agriculture, and natural risks analysis. Since few years, some high resolution DEMs have become open and freely available. Some of them are available on the Earth Engine data catalog, but only few are available at a very high resolution. And when it comes to analyse the human and financial consequences of natural disasters such as floods and sea level rise, this high resolution is needed.

Hence, the aim of this article is to detail the implementation steps of the high resolution France DEM (RGE ALTI® 5m — IGN) into Earth Engine. In the following, the procedure is applied on the southern part of Corsica (Department 2A) but can be easily transposed to the all country.

The structure of the article reads as follow:

• the dataset (RGE ALTI® 5m) is described,
• a procedure to automate the download of the DEM of multiple departments of France is developed from the FTP server,
• the DEM asc files are turned into geoTIFF files and merged into one,
• then, geoTIFF files are uploaded as new assets in Earth Engine.

The final result has been published on Awsome Google Earth Engine Community Datasets.

Run me first

The following librairies are required:

import os
import shutil

from ftplib import FTP

import py7zr

import pandas as pd
import numpy as np

import rasterio
from rasterio.merge import merge
from rasterio.plot import show

from osgeo import gdal, osr

import pprint
from tqdm.notebook import tqdm

import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
%matplotlib inline

# Prepare data folder
DATAPATH = 'data'
os.makedirs(DATAPATH, exist_ok=True)

Description of the RGE ALTI® 5m dataset

The RGE ALTI® 5m dataset describes the ground elevation of France with a spatial resolution of 5x5m. It is produced by the National Institute of Geographic and Forest Information (IGN — https://www.ign.fr/) which is a public administrative establishment placed under the joint authority of the Ministries in charge of ecology and forestry. The full dataset description is available here.

The accuracy of the DEM varies according to areas and issues: it is increased to 20 cm in floodplain areas or coastal areas. The RGE ALTI® is updated from surveys obtained by airborne LIDAR or by correlation of aerial images. The dataset can be downloaded at a 1x1m resolution or at a 5x5m resolution.

This dataset is published under the open licence Etalab.

Dataset accessibility — automate download and extraction

The RGE ALTI® 5m dataset is split in multiple 7-Zip files, and all individual files are available here. These files are stored on a FTP server ftp3.ign.fr with a public username/password.

Then, a procedure is given below to get the info and download files stored on the FTP associated with the dataset:

# Define requested info to access the FTP.
URL = "ftp3.ign.fr"
USERNAME = "RGE_ALTI_ext"
PASSWORD = "Thae5eerohsei8ve"

# Create a small helper class to manipulate our ftp connection
# without leaving the connection open between calls
class FTPHelper:
    def __init__(self, url, username, password):
        self.url = url
        self.username = username
        self.password = password
    
    def list_files(self):
        """
        Connects to ftp server and returns a list of file names stored on it.
        """
        with FTP(self.url) as ftp:
            ftp.login(self.username, self.password)
            filenames = ftp.nlst()
        return filenames


    def download_file(self, filename, local_folder, verbose=True):
        """
        Connects to ftp server and downloads remote `filename` to local `local_folder`.

        Disable verbose to avoid printing the progress bar.
        """
        # Create the directory to store the downloaded file.
        os.makedirs(local_folder, exist_ok = True)
        # Define the output filename.
        output_filepath = os.path.join(local_folder, filename)
        # Download the file (big files so ze show a progress bar)
        with FTP(self.url) as ftp:
            ftp.login(self.username, self.password)
            filesize = ftp.size(filename)
            with open(output_filepath, 'wb') as f:
                with tqdm(total=filesize,
                    unit='B', unit_scale=True, unit_divisor=1024,
                    disable=not verbose) as pbar:
                    pbar.set_description(f"Downloading {filename}")
                    def callback_(data):
                        l = len(data)
                        pbar.update(l)
                        f.write(data)
                    ftp.retrbinary('RETR ' + filename, callback_)
        print(f"{filename} is downloaded.")
        return output_filepath

# Instantiate helper
ftp_h = FTPHelper(URL, USERNAME, PASSWORD)

The FTPHelper.list_files() function is used to store all files info into a pandas dataframe.

# Create a pandas.df to store the list of available files.
df = pd.DataFrame(data=ftp_h.list_files(), columns=["name"])

# Rearrange the content of the pandas df.
df = pd.concat([df, df["name"].str.split("_", expand=True)], axis=1)
df = df.rename(columns = {
    0:"dataset", 1:"version", 2:"res", 3:"filetype", 4:"crs", 
    5:"dep", 6:"date"})

# Loc the dataset of interest RGE Alti 5m.
df = df.loc[(df.dataset == "RGEALTI") & (df.res == "5M")]

# Clean the date and set the name as index.
df['date'] = df['date'].apply(lambda x:x[:-3])
df = df.set_index("name")

df.head()

In the following we define a department of interest corresponding to the southern part of Corsica (France), and we print the associated file:

deps_of_interest = ['D02A']

# Loc lines of interest in the pandas.df.
dfi = df.loc[df.dep.isin(deps_of_interest)]
dfi

The file of interest is downloaded using our previously define function FTPHelper.download_file(). Please note that the server can be very slow.

filename = dfi.index[0]

# Define the name (global variable) of the temp directory to store zip files.
TEMP_ZIP = os.path.join(DATAPATH, "temp7z")

# Comment after running to no download again.
ftp_h.download_file(filename, TEMP_ZIP)

Then we need to extract the content of the downloaded file. A procedure is defined accordingly:

# Define the name (global variable) of the temp directory to store extracted files.
EXTRACTION_PATH = os.path.join(DATAPATH, "RGE_ALTI_5m")

def extract_rge5m(filename):
    """
    Extracts the zip DEM.
    """
    zipfile_path = os.path.join(TEMP_ZIP, filename)
    os.makedirs(EXTRACTION_PATH, exist_ok = True)
    with py7zr.SevenZipFile(zipfile_path, mode='r') as archive:
        archive.extractall(path=EXTRACTION_PATH) #WIN: check the path lenght in case of error FileNotFound
    print(f"{zipfile_path} is extracted to {EXTRACTION_PATH}.")

The function is applied to our filename of interest.

extract_rge5m(filename)

Now the unzip folder should look like this:

RGE_ALTI_5m
└───RGEALTI_2-0_5M_ASC_LAMB93-IGN78C_D02A_2020-04-16
    └───RGEALTI
        │   LISEZ-MOI.txt
        └───1_DONNEES_LIVRAISON_2021-06-00061
        │   │   RGEALTI_MNT_5M_ASC_LAMB93_IGN78C_D02A.md5
        │   └───RGEALTI_MNT_5M_ASC_LAMB93_IGN78C_D02A
        │       │   RGEALTI_FXX_1155_6135_MNT_LAMB93_IGN78C.asc
        │       │   RGEALTI_FXX_1155_6140_MNT_LAMB93_IGN78C.asc
        │       │   RGEALTI_FXX_1155_6145_MNT_LAMB93_IGN78C.asc
        │       │   ...
        └───2_METADONNEES_LIVRAISON_2021-06-00061
        └───3_SUPPLEMENTS_LIVRAISON_2021-06-00061

The folder/subfolder names depend on the acquisition period of the dataset, and on the department of interest. Hence, a procedure is given to go straight to the path where the asc files of the DEM are located. We store all asc files into a list. We finally print the 3 first asc paths of the list.

def get_path_asc_paths(filename):
    """
    This funciton returns a list of paths of asc files.
    """
    for root, dirs, files in os.walk(EXTRACTION_PATH, filename[:-3]):
        if "1_DONNEES_LIVRAISON" in root:
            asc_paths_list = sorted([os.path.join(root, name) for name in files if name.endswith(".asc")])

    return asc_paths_list

asc_paths_list = get_path_asc_paths(filename)
pprint.pprint(asc_paths_list[:3])

Let’s have a look at a DEM piece (.asc file) located in this folder. The 10th element is selected and the header of the file (6 first lines) is printed.

# Select one asc file in the DEM directory.
sample_path = asc_paths_list[10]

# Open and print the 10 first lines.
with open(sample_path) as file:
    content = file.readlines()
    header = content[:6]
    data = content[6:]

pprint.pprint(header)

As we can see, the header of the file gives us:

• ncols the number of columns,
• nrows the number of rows,
• xllcorner the x location (in epsg 2154) of the lower left corner of the DEM,
• yllcorner the y location (in epsg 2154) of the lower left corner of the DEM,
• cellsize the cell size in meters,
• NODATA_value the nodata value.

On the other hand, the rest of the content represents the elevation data. The 1000 first characters of the data content is printed as follow:

pprint.pprint(data[0][:1000])

We define a function to store these data into a dictionary. For further manipulations, it will be convenient to define a class associated with RGE pieces. Hence, we also define such a class in which the DEM data is loaded using the np.loadtxt() function.

def get_header_asc(filepath):
    """
    This function reads the header of an asc file and returns 
    the data into a dictionnary
    """
    file = open(filepath)
    content = file.readlines()[:6]
    content = [item.split() for item in content]
    return dict(content)

class RGEitem():
    """
    This class is used to handle RGE items.
    """
    def __init__(self, filepath):
        self.filename = os.path.basename(filepath)
        self.dir = os.path.dirname(filepath)
        self.data = np.loadtxt(filepath, skiprows=6)
        self.header = get_header_asc(filepath)
        self.ncols = int(self.header['ncols'])
        self.nrows = int(self.header['nrows'])
        self.xllc = float(self.header['xllcorner'])
        self.yllc = float(self.header['yllcorner'])
        self.res = float(self.header['cellsize'])
        self.zmin = float(self.data.min())
        self.zmax = float(self.data.max())
        self.novalue = -99999.

Let’s have a look at our sample item:

file = RGEitem(sample_path)

plt.rcParams["figure.figsize"] = (7,7)
plt.imshow(file.data, vmin=0, vmax=300)
plt.show()

print("The max elevation is", file.zmax,"m.")

The max elevation is 619.01 m.

Turn .asc files into geoTIFF files

Now that we see how to download and access these data, we can transform each .asc file into a geoTIFF using the osgeo library. We do this as follow:

def asc_to_tif(file, output_raster_dir, epsg):
    """
    Transforms an .asc file into a geoTIFF.

    Params:
    -------
    file: an RGEitem
    output_raster_dir (str): path to the directory where the tif will be saved.
    epsg (int): projection system.

    Returns:
    --------
    output_rasterpath (str): name of the output geoTIFF
    """
    xmin = file.xllc
    ymax = file.yllc + file.nrows * file.res 
    geotransform = (xmin, file.res, 0, ymax, 0, -file.res)

    output_rasterpath = os.path.join(output_raster_dir, file.filename[:-4] + ".tif")

    # Open the file
    output_raster = gdal.GetDriverByName('GTiff').Create(output_rasterpath, file.ncols, file.nrows, 1, gdal.GDT_Float32)
    # Specify the coordinates.  
    output_raster.SetGeoTransform(geotransform)
    # Establish the coordinate encoding.  
    srs = osr.SpatialReference()  
    # Specify the projection.               
    srs.ImportFromEPSG(epsg)                     
    # Export the coordinate system to the file.
    output_raster.SetProjection(srs.ExportToWkt())
    # Writes the array.   
    output_raster.GetRasterBand(1).WriteArray(file.data)
    # Set nodata value.  
    output_raster.GetRasterBand(1).SetNoDataValue(file.novalue) 
    output_raster.FlushCache()
    return output_rasterpath

We apply the function to get the geoTIFF associated with our sample and we open it using rasterio.

# Generate the tif.
SAMPLEPATH = os.path.join(DATAPATH, "sample")
os.makedirs(SAMPLEPATH, exist_ok = True)
tiffile = asc_to_tif(file, SAMPLEPATH, 2154)

with rasterio.open(tiffile) as img:
    show(img)
    plt.show()

Note that x and y axes are now in line with the coordinate reference system we are using (epsg:2154).

Merge geoTIFF files into a mosaic

We can now define a procedure to trnasform each .asc file into .tif and create a mosaic with all these files following the procedure described by Abdishakur.

def create_rge_mosaic(asc_paths_list, result_path, mosaic_name, crs):
    """
    Creates a mosaic associated to multiple asc files.

    Params:
    -------
    asc_paths_list (list): list of asc files paths.
    mosaic_path (str): path of the output mosaic.
    crs (int): coordinate reference system (ex. 2154 for EPSG:2154).
    """

    # Create tmp dir to save intermediate tifs.
    tmpdir = os.path.join(DATAPATH, 'local_tifs')
    os.makedirs(tmpdir, exist_ok = True)
    
    output_raster_paths_list = [asc_to_tif(RGEitem(ascpath), tmpdir, crs) for ascpath in asc_paths_list]
    
    # Safely load rasters and create the mosaic.
    from contextlib import ExitStack
    with ExitStack() as stack:
        raster_to_mosaic_list = [stack.enter_context(rasterio.open(path)) 
            for path in output_raster_paths_list
        ]
        mosaic, output = merge(raster_to_mosaic_list)
        output_meta = raster_to_mosaic_list[0].meta.copy()
    
    output_meta.update({
        "driver": "GTiff",
        "height": mosaic.shape[1],
        "width": mosaic.shape[2],
        "transform": output,
    })

    # Save the result.
    os.makedirs(result_path, exist_ok=True)
    mosaic_path = os.path.join(result_path, mosaic_name)
    with rasterio.open(mosaic_path, "w", **output_meta) as m:
        m.write(mosaic)
    
    # Purge the tmp dir.
    shutil.rmtree(tmpdir)

    return mosaic_path

The function is applied to create the mosaic of the southern part of Corsica (department 2A). The procedure an take few minutes and the final resulting file is about 1.4 GB:

result_path = os.path.join(DATAPATH, "result")
mosaic_name = "mymosaic2A_2022.tif"

mosaic_path = create_rge_mosaic(asc_paths_list, result_path, mosaic_name, 2154)

The result can be displayed as follows. On the following figure, the initial sample area is delineated with a red rectangle.

fig, ax = plt.subplots(1, figsize=(12, 12))

# Indicate the location of the sample tiff with a red rectangle.
x1, y1 = file.xllc, file.yllc
dx, dy = file.ncols*file.res, file.nrows*file.res
ax.add_patch(Rectangle((x1, y1), dx, dy, edgecolor='red', fill=False))

# PLot the mosaic of the 2A departement.
with rasterio.open(mosaic_path) as img:
    show(img, interpolation='none', ax=ax)

plt.title('Elevation of the southern part of Corsica (2A)')
plt.show()

Group mosaics and upload in Earth Engine as new assets

The procedure described in this notebook has been applied for all departments of France (Continental). Then, the resulting geoTIFFs have been uploaded into Google Earth Engine as new assets. The procedure to manage an assets in Earth Engine is described here.

The final result has been published on Awsome Google Earth Engine Community Datasets.

The Earth Engine snippet reads as follows:

var rge_alti5 = ee.Image("projects/sat-io/open-datasets/IGN_RGE_Alti_5m");

Sample Code: https://code.earthengine.google.com/b9054b6d5c4d7de1a4af39ff307e73d7

Acknowledgements

Thanks to Samapriya Roy for final merging of the dataset into Earth Engine, for publication on the AGCD, and for the production of the GIF used before intro.

Note

You can find this notebook and many others on my GitHub geoscience-notebooks repository or on my website.