pfg4
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
!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 |
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
# 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() |
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
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') |
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
# 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() |
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
# 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() |
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
# 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() |
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
# 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() |
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
# 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() |
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
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() |
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
# 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)) |
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
mymap.save("my_Python_for_Geosciences_map.html") |
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
# 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 |
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
dep["num"] = dep["code_insee"] | |
dep = dep.set_index("code_insee") | |
dep.at["69M", "num"] = "69" | |
dep.at["69D", "num"] = "69" |
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 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 |
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
myurl = getHubeauURL_piezo(mycode, mydate_i, mydate_f) | |
print(myurl) |
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
# We get the content of the webpage | |
r = requests.get(myurl) | |
# The server response: | |
print(r) |
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
# We organize the content into a json dictionnary | |
res = r.json() | |
# We print the result | |
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
res['data'][0] |
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 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() |
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=(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() |
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
# 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 |
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
llc = pnts[0][0] | |
trc = pnts[0][2] |
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 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 |
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
# 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() |
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
new_url = getHubeauURLstations("piezo", "2020-01-01", [llc, trc]) | |
stations = StationsToDF(new_url) | |
stations = stations.set_index("code_bss") | |
stations |
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
# 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)) |
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
# 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() |
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
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 |
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
sandre.loc[(sandre.code_parametre_sandre == 1301) | | |
(sandre.code_parametre_sandre == 1340)] |
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 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) |
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
q_url = getHubeauURLstations("quality", "2020-01-01", [llc, trc]) | |
q_stations = StationsToDF(q_url) | |
q_stations.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
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() |
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
# 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)) |
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
# 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 |
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
# 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() |
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 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 |
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
url_ex = getHubeauURLquality(code_ex, | |
date_initial_ex, | |
date_final_ex, | |
params_ex) | |
print(url_ex) |
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
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 |
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, 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() |
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 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 |
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
# 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() |
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
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() |
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
# 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') |
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
# 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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment