Skip to content

Instantly share code, notes, and snippets.

@guiattard
Last active February 21, 2022 14:13
Show Gist options
  • Save guiattard/afbb11ac613e4090f2b92d7f1514d047 to your computer and use it in GitHub Desktop.
Save guiattard/afbb11ac613e4090f2b92d7f1514d047 to your computer and use it in GitHub Desktop.
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