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