pfg6-recharge-ee
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) | |
) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Evaluate local precipitation conditions. | |
local_pr = pr.getRegion(poi, scale).getInfo() | |
pprint.pprint(local_pr[:5]) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
pr_df = ee_array_to_df(local_pr, ["precipitation"]) | |
pr_df.head(10) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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]) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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]) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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]) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import ee | |
# Trigger the authentication flow. | |
ee.Authenticate() | |
# Initialize the library. | |
ee.Initialize() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import pandas as pd | |
import matplotlib.pyplot as plt | |
import numpy as np | |
import folium | |
import pprint | |
import branca.colormap as cm |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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] |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 © <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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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()) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
zr = ee.Image(0.5) | |
p = ee.Image(0.5) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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" | |
) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Define the initial time (time0) according to the start of the collection. | |
time0 = meteo.first().get("system:time_start") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
initial_image = initial_rech.addBands( | |
ee.Image([initial_apwl, initial_st, initial_pr, initial_pet]) | |
) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
image_list = ee.List([initial_image]) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
arr = rech_coll.getRegion(poi, scale).getInfo() | |
rdf = ee_array_to_df(arr, ["pr", "pet", "apwl", "st", "rech"]).sort_index() | |
rdf.head(12) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment