pfg4

pfg4-install-import.py

!apt install gdal-bin python-gdal python3-gdal
!pip install fiona shapely pyproj
!apt install python3-rtree 
!pip install git+git://github.com/geopandas/geopandas.git
!pip install descartes 
!pip install git+https://github.com/python-visualization/folium

import pandas as pd
import geopandas as gpd
import zipfile
import requests
import numpy as np
from pyproj import Transformer
import os
import datetime
import matplotlib as mpl
import matplotlib.pyplot as plt

%matplotlib inline

pfg4-run1.py

# The url where we can fide the shapefile
url_dep = "http://osm13.openstreetmap.fr/~cquest/openfla/export/departements-20180101-shp.zip"

# The name of the zip file
file_dep = "departements-20180101-shp.zip"

# Command to donwload the file at the given url
r = requests.get(url_dep)

# Then we open the file
open(file_dep, 'wb').write(r.content)

# We extract the content of the .zip file
with zipfile.ZipFile(file_dep, 'r') as unzip: unzip.extractall("dep.shp")

# we finally read the shapefile and make some cleaning
dep = gpd.read_file("dep.shp")

# We remove the zipfile
os.remove(file_dep)

# we print the head of our geodataframe
dep.head()

pfg4-run10.py

dir_res = './RESULTS'
try:
# Create target Directory
    os.mkdir(dir_res)
except FileExistsError:
    pass
print('The folder is now created')

geo_bss.to_file('./RESULTS/all_boreholes.shp')
print('The file is saved')

pfg4-run11.py

# We create a new geodataframe with points where we know a groundwater table depth
geo_ww = geo_bss.loc[(geo_bss['point_eau'] == 'OUI') & (geo_bss['prof_eau_sol'].notna())]

# We only keep relevent features
geo_ww = geo_ww[['ID_BSS', 'code_bss', 'z_sol', 'z_bdalti','prof_investigation',
                 'prof_eau_sol', 'geometry']]
                 
# We create a new geodataframe with points where a geotehrmal activity is known
geoth_points = geo_bss.loc[geo_bss['categorie_geothermie'].notna() |
                    geo_bss['procede_geothermique'].notna() |
                    geo_bss['usage_geothermie'].notna()]

# We plot the result
fig, ax = plt.subplots(figsize=(10,10))

# We plot our deps of interest with the Pastel 1 colormap
dep.plot(ax = ax, column = 'num', cmap = 'Pastel1', alpha = 0.75)

# We add the location of borehole with water level in blue
geo_ww.plot(ax =ax, color = 'blue', markersize=1, label = 'Borehole with a groundwater table level')

# We add the location of geothermal installations in red
geoth_points.plot(ax =ax, color = 'red', markersize=5, label = 'Geothermal installation')

# We add some cities to our map
main_cities.plot(ax =ax, color = 'orange', marker = 'o', markersize=100, label = 'City')

plt.legend()
plt.show()

pfg4-run12.py

# we define the long./lat. of our location of interest
site_lon, site_lat = 4.8400, 45.7600 

site = pd.DataFrame({'Name': ['Lyon'],
                     'lat': [site_lat],
                     'lon': [site_lon]})

# We then build a geopandas dataframe to create our buffer zone:
buff_zone = gpd.GeoDataFrame(site, geometry = gpd.points_from_xy(site.lon, site.lat), 
                          crs = 'epsg:4326')

# to apply a buffer zone, we need to change our projeciton system to work in a metric system.
# In France, the common metric projection system is "epsg:2154"

buff_zone = buff_zone.to_crs('epsg:2154')

# we create our buffer zone of 10km:
buff_zone['geometry'] = buff_zone['geometry'].buffer(10000)

# We get back into epsg:4326 once our buffer zone is defined
buff4326 = buff_zone.to_crs('epsg:4326')
buff4326 = buff4326[['geometry']]

# We finally keep water wells informations located inside this studied area
geo_ww = gpd.sjoin(geo_ww, buff4326, op='within')

# We keep the columns we want:
geo_ww = geo_ww[['ID_BSS', 'code_bss', 'z_sol', 'z_bdalti','prof_investigation',
                 'prof_eau_sol', 'geometry']]

# We plot the result
fig, ax = plt.subplots(figsize=(10,10))
ax.set_xlim((min(geo_ww['geometry'].x), max(geo_ww['geometry'].x)))
ax.set_ylim((min(geo_ww['geometry'].y), max(geo_ww['geometry'].y)))

# We plot our deps of interest with the Pastel 1 colormap
dep.plot(ax = ax, column = 'num', cmap = 'Pastel1', alpha = 0.75)

# We add the location of borehole with water level in blue
geo_ww.plot(ax =ax, color = 'blue', markersize=1, label = 'Borehole with a groundwater table level')

# We add the location of geothermal installations in red
geoth_points.plot(ax =ax, color = 'red', markersize=5, label = 'Geothermal installation')

# We add our city of interest
ax.scatter(site_lon, site_lat, color = 'orange', marker = 'o', s=100, label = 'Point of interest')

plt.legend()
plt.show()

pfg4-run13.py

# First we convert all interesting fields to numeric values
for col in ['z_sol', 'z_bdalti', 'prof_investigation', 'prof_eau_sol']:
    geo_ww[col] = pd.to_numeric(geo_ww[col])

# Then we can describe our geodataframe
geo_ww.describe()

pfg4-run14.py

# We select rows inside the range we want
geo_ww = geo_ww[(geo_ww.z_sol > geo_ww.z_sol.quantile(0.01)) &
                (geo_ww.z_sol < geo_ww.z_sol.quantile(0.99))]

# Then we can describe our new geodataframe
geo_ww.describe()

pfg4-run15.py

# Calculation of the groundwater table elevation
geo_ww.loc[:, 'WT_elevation'] = geo_ww['z_bdalti'] - geo_ww['prof_eau_sol'] 

# We select wells with a relevant depth:
geo_ww = geo_ww[(geo_ww.prof_investigation >= 10) &
                (geo_ww.prof_investigation <= 30)]

# We remove nan values
geo_ww = geo_ww.dropna()

# We can now have a overview of different features of our dataset:
geo_ww[['prof_investigation', 'z_bdalti', 'WT_elevation']].hist()
#plt.tight_layout()
plt.show()

pfg4-run16.py

from mpl_toolkits.axes_grid1 import make_axes_locatable

# We plot the result
fig, ax = plt.subplots(figsize=(10,8))

ax.set_xlim((min(geo_ww['geometry'].x), max(geo_ww['geometry'].x)))
ax.set_ylim((min(geo_ww['geometry'].y), max(geo_ww['geometry'].y)))

# We plot our deps of interest with the Pastel 1 colormap
dep.plot(ax = ax, column = "num", cmap = "Pastel1", alpha = 0.75)

vmin = geo_ww['WT_elevation'].quantile(0.1)
vmax = geo_ww['WT_elevation'].quantile(0.9)

# We add the location of borehole with water level in blue
geo_ww.plot(ax = ax, cmap = 'viridis_r',
            column = 'WT_elevation',
            markersize=10, 
            vmin = vmin,
            vmax = vmax,
            legend = True,
            legend_kwds={'label': 'elevation [m a.s.l]'})

ax.set_title('Groundwater table elevation')
plt.show()

pfg4-run17.py

# We import the folium library
import folium

# We need to import branca.colormap to give pretty colors to our points according to groundwater
# table elevation
import branca.colormap as cm

# We build our map focusing on our site and specifying the zoom start
mymap = folium.Map(location=[site_lat, site_lon], 
                   zoom_start=12,
                   control_scale = True)

# Add the geological map of France
url = 'http://geoservices.brgm.fr/geologie'
layer1 = 'GEOLOGIE'
folium.WmsTileLayer(url, layer1, attr = 'BRGM', name = 'toto').add_to(mymap)
#folium.TileLayer('http://geoservices.brgm.fr/geologie/GEOLOGIE', attr = 'BRGM', name = 'toto').add_to(mymap)

# Add a feature group to add borehole with groundwater table elevation
fg_gwt = folium.FeatureGroup(name = 'Groundwater table elevation')
mymap.add_child(fg_gwt)

colormap = cm.LinearColormap(colors = ['orange', 'yellow', 'green', 'lightblue', 'blue'], 
                             index = [geo_ww["WT_elevation"].quantile(0.1),
                                    geo_ww["WT_elevation"].quantile(0.25),
                                    geo_ww["WT_elevation"].quantile(0.5),
                                    geo_ww["WT_elevation"].quantile(0.75),
                                    geo_ww["WT_elevation"].quantile(0.9)], 
                             vmin = geo_ww["WT_elevation"].quantile(0.1),
                             vmax = geo_ww["WT_elevation"].quantile(0.9))

# We add the caption of our colormap
colormap.caption = 'Groundwater table elevation [m asl]'

# We give an explicit description of longitude and latitude of our points:
geo_ww['lon'] = geo_ww['geometry'].x
geo_ww['lat'] = geo_ww['geometry'].y

#before interating over our dataset we reset its index
#geo_ww = geo_ww.reset_index()
# We add all points of our dataset
for i, h in zip(geo_ww.index.values, geo_ww.WT_elevation.values):
    dfi = folium.CircleMarker(
        location = [geo_ww.at[i,'lat'], geo_ww.at[i,'lon']],
        popup =str(round(h,2)) + ' m asl',
        radius= 7,
        fill=True,
        fill_opacity=0.7,
        color = colormap(h),
        fill_color = colormap(h))
    dfi.add_to(fg_gwt)

mymap.add_child(colormap)
mymap.add_child(folium.map.LayerControl(collapsed=True))

pfg4-run18.py

mymap.save("my_Python_for_Geosciences_map.html")

pfg4-run19.py

# Definition of an identifier:
mycode = "06987A0186/S"

# Definition of a period of interest:
mydate_i = "2015-01-01" #year/month/day
mydate_f = "2020-12-31" #year/month/day

pfg4-run2.py

dep["num"] = dep["code_insee"]
dep = dep.set_index("code_insee")
dep.at["69M", "num"] = "69"
dep.at["69D", "num"] = "69"

pfg4-run20.py

def getHubeauURL_piezo(code_bss, date_i, date_f):
    url_root = "https://hubeau.eaufrance.fr/api/v1/niveaux_nappes/chroniques?"
    url_code = "code_bss=" + code_bss.replace("/","%2F")
    url_date = "&date_debut_mesure=" + date_i + "&date_fin_mesure=" + date_f
    url_tail = "&size=20000&sort=asc"
    return url_root + url_code + url_date + url_tail

pfg4-run21.py

myurl = getHubeauURL_piezo(mycode, mydate_i, mydate_f)
print(myurl)

pfg4-run22.py

# We get the content of the webpage
r = requests.get(myurl)

# The server response:
print(r)

pfg4-run23.py

# We organize the content into a json dictionnary
res = r.json()

# We print the result
res

pfg4-run24.py

res['data'][0]

pfg4-run25.py

def piezoDico_to_pandas_df(res):
    # We convert the data into a dataframe
    df = pd.DataFrame.from_dict(res["data"])

    # We only keep relevant columns
    df = df[['date_mesure', 'niveau_nappe_eau']]

    # We rename columns
    df = df.rename(columns = {'date_mesure':'date',
                    'niveau_nappe_eau': 'hydraulic_head'})

    # We make sure that the hydraulic head is numeric
    df['hydraulic_head'] = pd.to_numeric(df['hydraulic_head'])

    # We convert the date as a datetime format and we define it as index
    df['date'] = pd.to_datetime(df['date'], format =  "%Y-%m-%d")
    df = df.set_index("date").to_period('d')
    return df

df = piezoDico_to_pandas_df(res)
df.head()

pfg4-run26.py

fig, ax = plt.subplots(figsize=(10,5))
plt.rcParams.update({'font.size': 12})

ax.set_title('Hydraulic head fluctuation at ' + mycode)

# We plot the mean hydraulic head over the period
df['hydraulic_head'].plot(ax = ax, 
              marker='o', 
              alpha=0.5,
              linestyle='--')

ax.set_ylabel('Hydraulic head (m)')
ax.set_xlabel('')

plt.legend()
plt.show()

pfg4-run27.py

# We define our rectangle of interest:
rect = gpd.GeoDataFrame(geometry = buff4326.envelope)

# We take put the coordinates into a list:
pnts = np.dstack(rect.at[0, 'geometry'].boundary.xy)
pnts

pfg4-run28.py

llc = pnts[0][0]
trc = pnts[0][2]

pfg4-run29.py

def getHubeauURLstations(param, date_research, bbox):
    """
    The aim of this function is to return the URL giving us available stations
    depending on:
    - param (str): the type of data we are looking for: param = "piezo" or "quality"
    - date_research (str): the date around when we want data ex. = "2020-01-01"
    - bbox: the area of research (list of two points coordinates). The fist element of the fist
    regroups coordinates of the lower left corner and the second element the top right corner.
    the function returns a URL
    """
    #param is "piezo" or "quality"
    if param == "piezo":
        url_root = "https://hubeau.eaufrance.fr/api/v1/niveaux_nappes/stations?"
    elif param == "quality":
        url_root = "https://hubeau.eaufrance.fr/api/v1/qualite_nappes/stations?"
    else:
        return print("error with param")

    llc = bbox[0]
    trc = bbox[1]
        
    llc_x, llc_y = str(round(llc[0],4)), str(round(llc[1],4))
    trc_x, trc_y = str(round(trc[0],4)), str(round(trc[1],4))
    
    url_bbox = "bbox=" + llc_x + "%2C" + llc_y + "%2C" + trc_x + "%2C" + trc_y
    url_date = "&date_recherche=" + date_research
    url_tail = "&format=json&size=20000"
    return url_root + url_bbox + url_date + url_tail

def StationsToDF(url):
    """
    The aim of this function is to organize the data of the resquested stations
    into a pandas df.
    The function takes an URL as input (str) and returns a pandas df
    """
    # We get the content of the webpage
    r = requests.get(url)
    # We organize the content into a json dictionnary
    res = r.json()
    # We convert the data into a dataframe
    df = pd.DataFrame.from_dict(res["data"])
    return df

pfg4-run3.py

# We select the numbers (as strings) of our department of interest
deps_of_interest = ['69', '01', '38']

# We reduce our geodataframe to our departement of interest
dep = dep.loc[dep["num"].isin(deps_of_interest)]
dep.plot()

plt.show()

pfg4-run30.py

new_url = getHubeauURLstations("piezo", "2020-01-01", [llc, trc])
stations = StationsToDF(new_url)

stations = stations.set_index("code_bss")
stations

pfg4-run31.py

# We import the folium library
import folium

# We build our map focusing on our site and specifying the zoom start
mymap = folium.Map(location=[site_lat, site_lon], 
                   zoom_start=11,
                   show=True, control_scale = True,
                   width = 400)

# Add the geological map of France
url = 'http://geoservices.brgm.fr/geologie'
layer1 = 'GEOLOGIE'
folium.WmsTileLayer(url, layer1, attr = 'BRGM', name = 'Geological map', show=True).add_to(mymap)

# Add a feature group to add borehole with active piezometric stations
fg_pz = folium.FeatureGroup(name = 'Active piezometric stations', show = True)
mymap.add_child(fg_pz)

# We add all points of our dataset
for code_bss in stations.index.values:
    pzi = folium.Marker(
        location = [stations.at[code_bss,'y'], stations.at[code_bss,'x']],
        popup = str(code_bss))
    pzi.add_to(fg_pz)

mymap.add_child(folium.map.LayerControl(collapsed=True))

pfg4-run32.py

# Definition of a period of interest:
mydate_i = "2015-01-01" #year/month/day
mydate_f = "2020-01-01" #year/month/day

fig, ax = plt.subplots(figsize=(7,5))
plt.rcParams.update({'font.size': 12})

for code in stations.index:
    myurl = getHubeauURL_piezo(code, mydate_i, mydate_f)
    r = requests.get(myurl)
    res = r.json()
    df = piezoDico_to_pandas_df(res)
    

    # We plot the mean hydraulic head over the period
    df['hydraulic_head'].plot(ax = ax, 
                  marker='o', 
                  alpha=0.5,
                  linestyle='--',
                  label = code)

ax.set_ylabel('Hydraulic head (m)')
ax.set_xlabel('')

plt.legend()
plt.show()

pfg4-run33.py

url_root= "http://www.sandre.eaufrance.fr/sites/default/files/document-sandre/"
file = "correspondance_par_sise_sandre_18092018.xls"
url = url_root + file

r = requests.get(url)

open(file, 'wb').write(r.content)
sandre = pd.ExcelFile(file).parse(0)
new_header = sandre.iloc[0]
sandre = sandre[1:]
sandre.columns = new_header
sandre = sandre[['Nom_famille_sise', 'nom_parametre_sise', 'code_parametre_sandre']]

sandre

pfg4-run34.py

sandre.loc[(sandre.code_parametre_sandre == 1301) |
           (sandre.code_parametre_sandre == 1340)]

pfg4-run35.py

def getSandreURL(code_param):
    """
    This function prints the SANDRE information associated to the code (input)
    """
    url_head = "http://www.sandre.eaufrance.fr/urn.php?urn=urn:sandre:donnees:PAR:FRA:code:"
    code_param = str(code_param)
    url_tail = ":::referentiel:2:html"
    url = url_head + code_param + url_tail
    print(url)
    
# Example:
getSandreURL(1301)

pfg4-run36.py

q_url = getHubeauURLstations("quality", "2020-01-01", [llc, trc])
q_stations = StationsToDF(q_url)
q_stations.head()

pfg4-run37.py

from datetime import datetime

q_stations['date_debut_mesure'] = pd.to_datetime(q_stations['date_debut_mesure'], format =  "%Y-%m-%d")
q_stations['date_fin_mesure'] = pd.to_datetime(q_stations['date_fin_mesure'], format =  "%Y-%m-%d")

datetime_i = datetime.strptime(mydate_i, "%Y-%m-%d")
datetime_f = datetime.strptime(mydate_f, "%Y-%m-%d")

q_stations = q_stations.loc[(q_stations['date_debut_mesure'] <= datetime_i) &
                            (q_stations['date_fin_mesure'] >= datetime_f)]

q_stations = q_stations.set_index("code_bss")
q_stations.head()

pfg4-run38.py

# We build our map focusing on our site and specifying the zoom start
mymap = folium.Map(location=[site_lat, site_lon], 
                   zoom_start=11,
                   show=True, control_scale = True,
                   width = 400)

# Add the geological map of France
url = 'http://geoservices.brgm.fr/geologie'
layer1 = 'GEOLOGIE'
folium.WmsTileLayer(url, layer1, attr = 'BRGM', name = 'Geological map', show=True).add_to(mymap)

# Add a feature group to add borehole with active quality stations
fg_quality = folium.FeatureGroup(name = 'Active groundwater quality stations', show = True)
mymap.add_child(fg_quality)

# We add all points of our quality dataset
for code_bss in q_stations.index.values:
    qi = folium.Marker(
        location = [q_stations.at[code_bss,'latitude'], q_stations.at[code_bss,'longitude']],
        popup = str(code_bss),
        icon=folium.Icon(color='red', icon='info-sign'))
    qi.add_to(fg_quality)
    
# Add a feature group to add borehole with active piezometric stations
fg_pz = folium.FeatureGroup(name = 'Active piezometric stations', show = True)
mymap.add_child(fg_pz)

# We add all points of our dataset
for code_bss in stations.index.values:
    pzi = folium.Marker(
        location = [stations.at[code_bss,'y'], stations.at[code_bss,'x']],
        popup = str(code_bss))
    pzi.add_to(fg_pz)

mymap.add_child(folium.map.LayerControl(collapsed=True))

pfg4-run39.py

# Definition of an identifier:
code_ex = "08036X1858/F2"

# Definition of a period of interest:
date_initial_ex = "2017-01-01" #year/month/day
date_final_ex = "2019-12-31" #year/month/day

params_ex = [1301, 1340] #our parameters in a list

pfg4-run4.py

# The url where we can fide the shapefile cities
url_cit = "https://simplemaps.com/static/data/country-cities/fr/fr.csv"

# Same procedure as seen before for departments
file_cit = "fr.csv"

r = requests.get(url_cit)
open(file_cit, 'wb').write(r.content)

cit = pd.read_csv('fr.csv', encoding='utf-8')

#We convert it into a geodataframe
geo_cit = gpd.GeoDataFrame(cit, geometry = gpd.points_from_xy(cit.lng, cit.lat),
                           crs = 'epsg:4326')

#We only keep cities inside our dep of interest:
geo_cit = gpd.sjoin(geo_cit, dep, op='within')

# We sort the associated table by population:
geo_cit = geo_cit.sort_values(by=['population_proper'], ascending = False).dropna()

# we print the head of our geodataframe
geo_cit.head()

pfg4-run40.py

def getHubeauURLquality(code, date_initial, date_final, params):
    """
    this function return the url needed to get the required parameters
    of the borehole code
    """
    url_head = "https://hubeau.eaufrance.fr/api/v1/qualite_nappes/analyses?"
    url_code = "bss_id=" + str(code).replace("/","%2F")
    url_params = "&code_param=" 
    for par in params:
        url_params += str(par) + "%2C%20"
        
    url_params = url_params[:-6]
    url_date = "&date_debut_prelevement=" + date_initial + "&date_fin_prelevement=" + date_final
    url_end = "&size=20000&sort=asc"
    
    url = url_head + url_code + url_params + url_date + url_end

    return url

pfg4-run41.py

url_ex = getHubeauURLquality(code_ex, 
                             date_initial_ex, 
                             date_final_ex, 
                             params_ex)

print(url_ex)

pfg4-run42.py

r = requests.get(url_ex)
res = r.json()
df = pd.DataFrame.from_dict(res["data"])

df = df.rename(columns={"date_debut_prelevement":"date",
                        "resultat": "value"})

df = df[['date', 'code_param', 'nom_param', 'value',
         'limite_quantification','limite_detection', 'incertitude_analytique']]

# We explicitly define the date in a datetime format
df['date'] = pd.to_datetime(df['date'], format =  "%Y-%m-%d")
df = df.set_index("date")
df

pfg4-run43.py

fig, axs = plt.subplots(1, 2, figsize=(15,5))
plt.rcParams.update({'font.size': 12})
ax = axs[0]
ax.set_title('Temperature fluctuation at ' + mycode)

# We plot the temperature fluctuation over the period
df.loc[df.code_param == 1301]['value'].plot(ax = ax, 
              marker='o', 
              alpha=0.5)

ax.set_ylabel('Temperature [°C]')
ax.set_xlabel('')

ax = axs[1]
ax.set_title('Nitrates concentration fluctuation at ' + mycode)

# We plot the nitrates fluctuation over the period
df.loc[df.code_param == 1340]['value'].plot(ax = ax, 
              marker='o', 
              alpha=0.5)

ax.set_ylabel('Concentration [mg/L]')
ax.set_xlabel('')

plt.show()

pfg4-run5.py

def download_bss_dep(num):
    # We use the URL where we can find all BSS zipfile sorted by department
    url_root = "http://infoterre.brgm.fr/telechargements/ExportsPublicsBSS/"
    file = "bss_export_" + num + ".zip"
    url = url_root + file
    
    # It will be the final path of our file
    final_path = "bss_export_" + num +".csv"
    
    if os.path.isfile(final_path):
        # If we have already download this BSS dep, no need to download it again
        print(file +" already exists. Remove it if an update is necessary.")
        bss = pd.read_csv(final_path, sep=';', decimal=",", low_memory=False)
    else:
        # We download the data and we reproduce the previous process
        print(file +" does not exist. Download in progress.")
        r = requests.get(url)
        open(file, 'wb').write(r.content)
        with zipfile.ZipFile(file, 'r') as unzip: unzip.extractall()
        print(file +" is now downloaded.")

        # We remove the zip
        os.remove(file)

        # We build the dataframe
        bss = pd.read_csv(final_path, sep=';', decimal=",", low_memory=False)
    return bss

pfg4-run6.py

# We create an empty dataframe
bss_dataset = pd.DataFrame()

# Our loop where we apply our function
for num_dep in deps_of_interest:
    bss = download_bss_dep(num_dep)
    frames = [bss_dataset,bss]
    bss_dataset = pd.concat(frames)
    print('Departement', num_dep, 'is done.')
    
bss_dataset.head()

pfg4-run7.py

items = ['ID_BSS', 'indice', 'designation', 'x_ref06', 'y_ref06','lex_projection_ref06','z_sol', 'z_bdalti', 'point_eau',
       'prof_investigation', 'prof_eau_sol', 'lex_usage', 'procede_geothermique',
       'categorie_geothermie', 'usage_geothermie']

# We reduce our dataframe to relevant items:
bss_dataset = bss_dataset[[col for col in items]]

# We need to recreate the "old" ID of boreholes to access some data (later)
bss_dataset['code_bss'] = bss_dataset['indice'] + '/' + bss_dataset['designation']
bss_dataset.head()

pfg4-run8.py

# We build a epsg df to associate projection names and epsg codes:
epsg_df = pd.DataFrame(columns = ['name', 'code_epsg'])
epsg_df = epsg_df.set_index('name')
epsg_df.at['Lambert-93', 'code_epsg'] = '2154'
epsg_df.at['Antilles-84-RRAF91', 'code_epsg'] = '32620'
epsg_df.at['Guyane-95', 'code_epsg'] = '2972'
epsg_df.at['Réunion-92', 'code_epsg'] = '2975'
epsg_df.at['Mayotte-2004', 'code_epsg'] = '4471'

# We create a longitude and latitude column to our dataset:
bss_dataset['long'] = ''
bss_dataset['lat'] = ''

# We make a loop through different projections to convert all in long/lat
for name in epsg_df.index:
    if name in set(bss_dataset['lex_projection_ref06'].dropna()):
        code = epsg_df.at[name, 'code_epsg']
        
        # This is the input projection name
        inProj = 'epsg:' + code
        
        # This is the output projection name
        outProj = 'epsg:4326'

        # We get the X and Y data
        x = np.asanyarray(bss_dataset.loc[bss_dataset['lex_projection_ref06'] == name]['x_ref06'])
        y = np.asanyarray(bss_dataset.loc[bss_dataset['lex_projection_ref06'] == name]['y_ref06'])
        
        # We define the transformer
        transformer = Transformer.from_crs(inProj, outProj)
        
        # We apply the transformer
        latitude ,longitude = transformer.transform(x, y)
        
        # We add the long/lat info to the dataframe
        bss_dataset.loc[bss_dataset['lex_projection_ref06'] == name, 'long'] = longitude
        bss_dataset.loc[bss_dataset['lex_projection_ref06'] == name, 'lat'] = latitude
    else:
        continue

# we remove lines where longitudes/latitudes are empty:
bss_dataset = bss_dataset.loc[(bss_dataset['long'] != '') &
                              (bss_dataset['lat'] != '')]

# Finally we can convert our dataframe into a geodataframe:
geo_bss = gpd.GeoDataFrame(bss_dataset, 
                           geometry=gpd.points_from_xy(bss_dataset.long, bss_dataset.lat),
                           crs = 'epsg:4326')

pfg4-run9.py

# We make some cleaning by removing points which are outside our departments of interest
geo_bss = gpd.sjoin(geo_bss, dep, op='within')

fig, ax = plt.subplots(figsize=(10,10))

# We plot our deps of interest with the Pastel 1 colormap
dep.plot(ax = ax, column = 'num', cmap = 'Pastel1', alpha = 0.75)

# We add the location of borehole in black
geo_bss.plot(ax =ax, color = 'black', markersize=1, label = 'Borehole', alpha = 0.25)

# We add 5 main cities
main_cities = geo_cit.head(5)
main_cities.plot(ax =ax, color = 'orange', marker = 'o', markersize=100, label = 'City')

# We add some cities to our map
for idx, row in main_cities.iterrows():
    txt = plt.annotate(s=row['city'], 
                 xy=[row['geometry'].x +0.025, row['geometry'].y +0.025],
                 horizontalalignment='left',
                 color = 'black',
                 backgroundcolor='w',
                 fontsize = 12)
    
    txt.set_bbox(dict(facecolor='w', alpha=0.5, edgecolor='w'))

plt.legend()
plt.show()