Introduction

In France, main geological/geoscientific data are managed by the French Geological Survey and several platforms provide a public and free access to many datasets such as:

  • geological maps,
  • boreholes, piezometers and wells locations, geological reports and cross sections (InfoTerre). In France, this database is known as the BSS (Base de donnée du Sous-Sol) which includes more than 700 000 boreholes descriptions,
  • piezometric time series and physico-chemical analysis from the national monitoring network (ADES).

These data dedicated to the (hydro-)geological knowledge are used by engineers and scientists to prevent natural risks and disasters (e.g. landslides or settlements), protect groundwater ressources and secure water supply for example. Particularly, these data can be used in many fields such as hydrogeology, geothermal energy and geotechnical applications.

These datasets can be explored and download manually by visiting the main platform of the french geological survey InfoTerre, but when it comes to export a large amount of data, at a regional scale, this task becomes tedious and time consuming. To face this problem and improve the dissemination of geoscientific data to the public and to scientists, the French geological survey recently developed many geo-webservices in accordance with European open-data requirements. These geo-webservices can help to automate geological data manipulation, analysis and mapping.

Consequently, the aim of this article is to present some of these services. Particularly, we are going to:

  • download and analyze boreholes information at the scale of the department of our choice: water supply wells and geothermal installations will be identified.
  • map all this hydrogeological information on an interactive folium map and add the geological tiles of France.
  • play with two Hubeau APIs to describe piezometric and physico-chemical characteristics of some groundwater wells: we will learn how to access and analyze hydraulic heads and groundwater quality fluctuations on the groundwater wells we want.

Run me first

Install/import usefull libraries

Before starting we need to install/import some libraries:

Administrative units of France

To download information at the scale of a department, we need some a geographical dataset with administrative units of France. We can use the openstreetmap shapefile of departments of France. We download this shapefile and put it in a local directory as follow:

code_insee nom nuts3 wikipedia surf_km2 geometry
0 974 La Réunion FR940 fr:La Réunion 2505 MULTIPOLYGON (((55.21643 -21.03904, 55.21652 -…
1 11 Aude FR811 fr:Aude (département) 6343 POLYGON ((1.68872 43.27368, 1.69001 43.27423, …
2 43 Haute-Loire FR723 fr:Haute-Loire 5003 POLYGON ((3.08206 45.28988, 3.08209 45.29031, …
3 13 Bouches-du-Rhône FR823 fr:Bouches-du-Rhône 5247 MULTIPOLYGON (((4.23014 43.46047, 4.23025 43.4…
4 47 Lot-et-Garonne FR614 fr:Lot-et-Garonne 5385 POLYGON ((-0.14058 44.22648, -0.12931 44.23218…

All department numbers are indicated in the column code_insee. We should be careful about this number because the Rhône department is identified as 69D to make the difference with the metropolitan area of Lyon which is identified as 69M. Here, we do not want to distinguish them so we create a new column associating the same number 69.

Selection of some departments of interest

To see the result, we can select some departments of interest and preview their shape. In the following, we will work with this selection of departments numbers to get some (hydro-)geological information. I decided to work with the department of the Ain (01), the Rhône (69) and of the Isère (38).

We also download a dataset of main cities of France (more than 100’000 inhabitants) and we keep the cities inside our depatment of interest following a similar procedure.

city lat lng country iso2 admin_name capital population population_proper geometry index_right
8 Lyon 45.76 4.84 France FR Auvergne-Rhône-Alpes admin 516092 516092 POINT (4.84000 45.76000)
5 Grenoble 45.1715 5.7224 France FR Auvergne-Rhône-Alpes minor 687985 158454 POINT (5.72240 45.17150)
165 Bourg-en-Bresse 46.2056 5.2289 France FR Auvergne-Rhône-Alpes minor 41527 41527 POINT (5.22890 46.20560)
256 Vienne 45.5242 4.8781 France FR Auvergne-Rhône-Alpes minor 29306 29306 POINT (4.87810 45.52420)

 

Extraction and analysis of boreholes information (BSS data)

Get our data into a clean (geo-)dataframe

Now, we want to get some information about the geological units of these departments. In the following, we are going to define a function that take the department number as input and which give a dataframe of boreholes information as output. Then, we will be able to identify we groundwater wells or geothermal installations are located and determine the groundwater table depth.

Now we can make a loop and concatenate the data of our departments of interest:

ID_BSS indice designation x_ref06 y_ref06 lex_projection_ref06 code_bss
0 BSS001PYDG 06246X0004 CPT 810400 6575119 Lambert-93 06246X0004/CPT
1 BSS001PYDH 06246X0005 CPT 810316 6574955 Lambert-93 06246X0005/CPT
2 BSS001PYDJ 06246X0006 CPT 810209 6575025 Lambert-93 06246X0006/CPT
3 BSS001PYER 06248X0003 S 827507 6575474 Lambert-93 06248X0003/S
4 BSS001PYES 06248X0004 HY 828134 6576416 Lambert-93 06248X0004/HY

 

As you can see, there are 71 columns. In the following, we reduce our focus to few relevent attributs (but please add more if you want to):

  • ID_BSS, indice and designation which are national codes of boreholes. These codes are requested when making queries to different webservices (we’ll see that later),
  • point_eau is a bolean indicating if there is some groundwater or not,
  • x_ref06 and y_ref06 are indicating the location in a geographical projection given by lex_projection_ref06,
  • z_sol indicates the ground elevation mesured and z_bdalti the ground elevation given by the IGN at the same location (can be useful to get both to identify some errors),
  • prof_investigation indicates the total depth of the borehole,
  • prof_eau_sol indicates the groundwater table depth if any exists,
  • lex_usage refers to the usage of the borehole/water well,
  • procede_geothermique, categorie_geothermie, and usage_geothermie are indicating if the point is actually a geothermal installation.

Your content goes here. Edit or remove this text inline or in the module Content settings. You can also style every aspect of this content in the module Design settings and even apply custom CSS to this text in the module Advanced settings.

Because the projection is not the same for all departments of France (because of overseas departments), the epsg number of different projections must be known. They are listed in a epsg dataframe below and finally, a procedure is made to convert all location in longitudes/latitudes (epsg:4326). The following step is required in case you want to work on overseas departments. In case you won’t, the following procedure still works to get longitudes/latitudes from Lambert 93 coordinates.

In practice, there are some errors in the reporting of points coordinates so to remove these errors, we make a spatial join of our geodataframe with the spatial extent of our department of interest: we only keep points inside our departments. Finally, we plot the location of boreholes.

daytime-land-surface-temperature-france

We can of course save our geodataframe as a shapefile, in a RESULTS folder as follow:

Identifying groundwater wells and geothermal installations

Of course we can decide to distinguish the kind of borehole we want to keep. For example, let’s select (1) boreholes where a groundwater table depth is reported and (2) geothermal installations.

In the first case, we keep points where the attribute point_eau is OUI (meaning that there is some groundwater right here) AND where prof_eau_sol is not a NaN value (meaning that we should know the groundwater table depth). In the second case, we keep points where categorie_geothermie or procede_geothermique or usage_geothermie are not a NaN value.

ground-elevation-france

Describing groundwater table elevation

Now that we know where groundwater bodies are located, we might want to describe groundwater levels and maybe describe groundwater flow directions. To do so, we refine a location of interest by choosing the coordinate of a point, we also define a relevant buffer zone around this point to reduce our dataset to this area of interest. In our example, we select the city of Lyon, and a buffer zone of 10km.

Please note that groundwater table elevation fluctuates over time. In the following, the calculated value represents the water table elevation observed when the borehole was drilled. Also, it should be noted that the calculated values are not synchronous which can alter the local representativity regarding the interpretation groundwater flow directions and flow systems.

Please note that we could instead import the polygon shapefile of our choice and select boreholes inside our polygon.

We describe data inside our geodataframe:

  z_sol z_bdalti prof_investigation prof_eau_sol
count 2080 1860 2008 2092
mean 174.73377 181.92688 19.682102 8.022075
std 86.678638 30.315466 13.985364 6.130222
min -999 155 2.3 0.2
25% 166.7975 167 12 4.2575
50% 169.92 170 17.5 6.1
75% 182 184 22.5 10
max 370 379 250.25 85

 

The result shows some errors regarding the ground elevation (-999 value). Consequently, let’s remove data outside first and last percentiles:

  z_sol z_bdalti prof_investigation prof_eau_sol
count 2036 1826 1953 2036
mean 178.90788 180.04272 19.775622 8.048355
std 23.567719 24.353793 14.11073 6.157885
min 157.044 155 2.3 0.2
25% 166.9 167 12 4.3
50% 169.945 170 17.5 6.1
75% 182 183 23 10
max 329 330 250.25 85

Now, we define the groundwater table elevation as the ground elevation minus the groundwater table depth. Because we have two different possibility to choose the ground elevation from the IGN z_bdalti. In addition, we can consider that in our area of interest we focus on the shallow aquifer where groundwater is flowing in the modern alluvial deposit. Consequently, let’s focus on an investigation depth ranging from 10m to 30m:

Then, groundwater table elevation values can be displayed as follow:

Interactive mapping with folium

The information above can easily be projected on an interactive map using the folium library. Additionally, it is interesting to do so because it can help us to understand the relation with geological units.

In fact, since a few years, geological tiles are available across France and can be projected and visualized using folium. All WMF/WFS services provided by the BRGM are listed in a table at this page. It is indicated that the geological map can be visualized using the following URL http://geoservices.brgm.fr/geologie, and the following layer GEOLOGIE. We directly use these inputs to define our WMS tile layer (see below in the script).

Then we can build our map using folium:

Of couse we can save this map in html and open it later with any navigator:

Access to hydraulic head monitoring data

Since a couple of years, several APIs dedicated to water management have been developed on a national platform called Hubeau. Some of these APIs are dedicated to hydrobiology and surface water monitoring, and others to quantitative and qualitative description of groundwater.

The principle is quite easy: you provide the API of your choice with (1) one or several borehole identifiers, (2) a period of interest, (3) some parameters of interest (piezometric level, temperature, etc.) and you get back a dictionary with all information.

First the procedure to describe the hydraulic head fluctuation at a given borehole will be described. Secondly, we will focus on how to find all recording stations around a area of interest, for a given period of time.

Getting hydraulic head fluctuations

Let’s see how the groundwater piezometric API works. We first need to define one borehole identifier code_bss, and a period of interest:

To query the API, we need to build an URL which depends on previous parameters. I give the function below providing you the appropriate URL:

We can now call this function as follow:

It returns the following URL you can visit. We can now make our query using the resquests.get function:

The code ‘200’ means that our query is successful (the code is valid and there are some data available). We now organize the result (using the json method) and print the content of the webpage using the following procedure:

The resulting content is a dictionary organized according to the date of measurements. We now need to define a function to extract the information we want in this dictionary. The first key of this dictionary gives us a count of all measurements over our period of interest. All piezometric values are stored in a data sub-dictionary. Let’s see on the first element of the data dictionary:

{‘code_bss’: ‘06987A0186/S’,
‘code_continuite’: ‘2’,
‘code_nature_mesure’: None,
‘code_producteur’: ‘196’,
‘date_mesure’: ‘2015-01-01’,
‘mode_obtention’: ‘Valeur mesurée’,
‘niveau_nappe_eau’: 163.09,
‘nom_continuite’: ‘Point lié au point précédent’,
‘nom_nature_mesure’: None,
‘nom_producteur’: ‘Service Géologique Régional Rhône-Alpes (196)’,
‘profondeur_nappe’: 5.41,
‘qualification’: ‘Correcte’,
‘statut’: ‘Donnée contrôlée niveau 2’,
‘timestamp_mesure’: 1420117200000,
‘urn_bss’: ‘
http://services.ades.eaufrance.fr/pointeau/06987A0186/S‘}

We have here a new dictionary with the ID of our borehole, the date of the measurement (date_mesure), the associated timestamp (timestamp_mesure), the hydraulic head measured (niveau_nappe_eau) and other information.

We can make a function to convert this dictionary into a dataframe keeping only relevant features as follow:

date hydraulic_head
01/01/2015 163.09
02/01/2015 163.08
03/01/2015 163.07
04/01/2015 163.06
05/01/2015 163.19

We can finally plot the hydraulic head fluctuation:

Finding piezometric monitoring stations around a location of interest

A procedure exists to find all stations with hydraulic head records around a location of interest. What you need is to define the longitude and latitude coordinates of a rectangle area. Let’s take the smallest rectangle including our previous circular buffer zone:

To make our query, we only need the lower left and top right corners coordinates:

We then need to make a new function to produce the URL and another one to organize the result in a pandas dataframe:

 

We run this function to sew how many piezometric stations are working/recording at a given date: the 1st of january 2020.

code_bss date_debut_mesure date_fin_mesure nom_commune x y bss_id
06988C0281/F 10/04/2009 17/08/2020 Chassieu 4.969186 45.742375 BSS001TPVW
07223X0130/P 18/05/2005 03/03/2020 Saint-Priest 4.908662 45.701409 BSS001URVP
06987X0272/P 14/09/2005 09/03/2020 Vaulx-en-Velin 4.937039 45.807487 BSS001TPEW
06987X0274/P 14/09/2005 09/03/2020 Vaulx-en-Velin 4.910798 45.796308 BSS001TPEY
06988X0217/P 23/08/2005 09/03/2020 Décines-Charpieu 4.952973 45.776655 BSS001TQGK
06987A0186/S 16/04/1968 07/08/2020 Villeurbanne 4.864856 45.779492 BSS001TMCR
06987J0105/PZ 16/05/2005 03/03/2020 Villeurbanne 4.899866 45.782911 BSS001TNEM

We can of course display the location of these active stations on an interactive folium map:

Getting hydraulic head fluctuations

Now, we can plot the piezometric levels of these stations over the similar period of interest that we defined earlier.

Access to physico-chemical parameters

Similarly, we can make some queries to observe the fluctuation of some physico-chemical parameters such as groundwater temperature, nitrates concentrations, pesticides… many other parameters are available.

Exploration of accessible physico-chemical parameters

Since 1992, the French producers of public water data have engaged a consistency process for their data in the framework of the French water information system. Consequently, all groundwater quality parameters are associated to a number in a referential SANDRE.

This referential can be download below. We can see that 1845 parameters are available.

 

Nom_famille_sise nom_parametre_sise code_parametre_sandre
1 COMPOSES ORGANOHALOGENES VOLATILS Tétrachloroéthane-1,1,1,2 1270
2 COMPOSES ORGANOHALOGENES VOLATILS 1,1,1,2 Tétrachloropropane 2704
3 COMPOSES ORGANOHALOGENES VOLATILS 1,1,1,3 Tétrachloropropane 2705
4 COMPOSES ORGANOHALOGENES VOLATILS Trichloroéthane-1,1,1 1284
5 COMPOSES ORGANOHALOGENES VOLATILS Tétrachloroéthane-1,1,2,2 1271
1841 PHYTOPLANCTONS % de colo de ulothricophyc. subst 6431
1842 PHYTOPLANCTONS Colonies de zygophyc. substrat 6432
1843 PHYTOPLANCTONS % de colo de zygophyc. substr 6432
1844 PARAMETRES MICROBIOLOGIQUES Entérocoques / 100 mL (qPCR) 6455
1845 PHYTOPLANCTONS Radiocystis sp (cellules) substrat 7299

We could explore all accessible parameters, but in this tutorial, we focus on two of them:

  • temperature with the code ‘1301’,
  • nitrates concentration with the code ‘1340’.

It can be checked by printing the appropriate lines of the dataframe:

  Nom_famille_sise nom_parametre_sise code_parametre_sandre
1255 PARAMETRES AZOTES ET PHOSPHORES Nitrates (en NO3) 1340
1645 CONTEXTE ENVIRONNEMENTAL Température de l’eau 1301

 

The information about these parameters (units, etc.) are described in the following webpages (in french):

You can easly access and get the information about the parameter you want by using the following function depending on the parameter code:

Finding quality monitoring stations around a location of interest

First, we find some quality monitoring stations around our location of interest by using the function defined earlier:

  bss_id code_bss date_debut_mesure date_fin_mesure longitude latitude altitude nom_commune
0 BSS001TLEA 06986T0059/S7 05/02/2004 05/02/2004 4.840702 45.731292 168 Lyon 7e Arrondissement
1 BSS001TPEA 06987X0252/PZ4 14/10/2002 02/06/2004 4.875022 45.762547 168 Villeurbanne
2 BSS001TPFA 06987X0276/PZ2 31/10/2003 04/10/2004 4.869878 45.750762 170 Lyon
3 BSS001URFA 07222X0490/PN4 27/11/2000 12/12/2007 4.845675 45.692831 164 Saint-Fons
4 BSS001URGA 07222X0514/PZ1 None None 4.845893 45.706085 164 Saint-Fons

 

We see that there are 471 stations. However, the dates are not necessarily matching our query (it is probably an issue that will be solved in a future version of the API). To refine our selection on our period of interest we convert all date features in datetimes format and compare with our datetimes of interest:

We have finally 7 active stations and we can add them (red color) to our interactive map:

Getting groundwater temperature and nitrates concentration fluctuations
We first need to define one borehole identifier code_bss, a period of interest and the code of parameters we want to explore (temperature and nitrates in our example):

To query the API, we need to build an URL which depends on previous parameters. I give the function below providing you the appropriate URL depending on the period and the parameters we want:

We can now call this function as follow:

Now, we get and organize the content of this webpage into a pandas dataframe. Please note that a very large number of features are available. You can explore these features by visiting the API website or by removing the last command in the following cell.

Here we only keep following features:

  • the date of the measurement,
  • the code of the parameter measured,
  • the name of the parameter measured (in french),
  • the resulting value,
  • the quantification limit,
  • the detection limit,
  • the analytical uncertainty

and we finally print our dataset:

date code_param nom_param value limite_quantification limite_detection incertitude_analytique
11/05/2017 1340 Nitrates 19 0.1 None None
07/12/2017 1340 Nitrates 19 0.1 None None
07/12/2017 1301 Température de l’Eau 16.7 NaN None None
05/07/2018 1340 Nitrates 17.5 0.1 None None
05/07/2018 1301 Température de l’Eau 17 NaN None None
16/10/2018 1301 Température de l’Eau 17 NaN None None
16/10/2018 1340 Nitrates 17.9 0.1 None None
21/05/2019 1340 Nitrates 20 0.25 None None
21/05/2019 1301 Température de l’Eau 17 NaN None None
16/10/2019 1301 Température de l’Eau 17.2 NaN None None
16/10/2019 1340 Nitrates 18 0.25 None None

 

We can now plot the results:

Discussion

This notebook intends to introduce some APIs and geo-webservices dedicated to the (hydro-)geological knowledge of France. This notebook has been reviewed and successfully returns several parameter descriptions around an area of interest. Please note that all input data have been carefully checked, and be aware that an alteration of input data might raise errors and/or inconsistencies. Particularly, some site-specific choices have been made to illustrate different possibilities and must not be generalized to other locations. Of course, readers are invited to play with all input parameters, and focus on other locations, but following advises and comments are given:

  • the area of interest should be carefully determined and must be of a reasonable size to limit the execution time of different procedures. Particularly, launching a data extraction of the all France is misadvised,
  • the closeness and depth of boreholes must be taken carefully when comparing some measurements: un-correlated fluctuations can be observed when sensors are note located into the same groundwater body or geological unit, even if boreholes are close.
  • there are many more available API functionalities not discussed in this notebook. Please visit Hubeau and PIC’EAU to find out more.
  • to work at reginal scales on several stations at a time, there are more efficient methods to extract and analyse data than ones presented in this notebook.

Play with the script…

This post is available as a notebook and can be run on Goolge Colab.

Please note that all my notebooks are available on my GitHub repository PythonForGeosciences.