130 lines
4.9 KiB
Python
130 lines
4.9 KiB
Python
import pandas as pd
|
|
import geopandas as gpd
|
|
import os
|
|
import pyproj
|
|
import matplotlib.pyplot as plt
|
|
from shapely.geometry import Point
|
|
|
|
def project_shp(L, wgs84_to_l93, l93_to_wgs84):
|
|
i = 0
|
|
while i < len(df):
|
|
if not df.loc[i, "has_gps"]:
|
|
start = i - 1
|
|
while i < len(df) and not df.loc[i, "has_gps"]:
|
|
i += 1
|
|
end = i
|
|
|
|
true_end = None
|
|
for j in range(end, min(end + 2000, len(df))):
|
|
if df.loc[j, "has_gps"] and df.loc[j, "Distance_entre_points_m"] > 0:
|
|
true_end = j
|
|
break
|
|
|
|
if true_end:
|
|
x_l93, y_l93 = wgs84_to_l93.transform(df.loc[start, "Longitude_WGS"],df.loc[start, "Latitude_WGS"])
|
|
P_start = Point(x_l93, y_l93, )
|
|
dist_start = L.project(P_start)
|
|
|
|
x_l93, y_l93 = wgs84_to_l93.transform(df.loc[true_end, "Longitude_WGS"],df.loc[true_end, "Latitude_WGS"])
|
|
P_end = Point(x_l93, y_l93)
|
|
dist_end = L.project(P_end)
|
|
|
|
n = true_end - start - 1
|
|
if n > 0:
|
|
step = (dist_end - dist_start) / (n + 1)
|
|
idx = start
|
|
# dist_temp_cumul = dist_start
|
|
for k in range(1, n + 1):
|
|
# dist_temp_cumul += step
|
|
# df.at[idx, "Distance_cumulee_corrigee"] = dist_temp_cumul
|
|
|
|
try:
|
|
x_l93, y_l93 = wgs84_to_l93.transform(df.loc[idx, "Longitude_WGS"],df.loc[idx, "Latitude_WGS"])
|
|
|
|
point = Point(x_l93, y_l93)
|
|
actual_distance = L.project(point)
|
|
|
|
projected_point = L.interpolate(actual_distance + step)
|
|
|
|
x_wgs, y_wgs = l93_to_wgs84.transform(projected_point.x, projected_point.y)
|
|
|
|
df.at[idx + 1, "Latitude_WGS"] = y_wgs
|
|
df.at[idx + 1, "Longitude_WGS"] = x_wgs
|
|
|
|
except:
|
|
continue
|
|
idx += 1
|
|
i += 1
|
|
|
|
|
|
def calcul_distance_cumulee(L, wgs84_to_l93, l93_to_wgs84):
|
|
|
|
x_l93, y_l93 = wgs84_to_l93.transform(df.loc[0, "Longitude_WGS"],df.loc[0, "Latitude_WGS"])
|
|
df.at[0, "Longitude_Lambert"] = x_l93
|
|
df.at[0, "Latitude_Lambert"] = y_l93
|
|
|
|
df.at[0, "Distance_Cele"] = L.project(Point(x_l93, y_l93))
|
|
df.at[0, "Distance_Cumulee"] = 0.0
|
|
|
|
for i in range(1, len(df)):
|
|
x_l93, y_l93 = wgs84_to_l93.transform(df.loc[i, "Longitude_WGS"],df.loc[i, "Latitude_WGS"])
|
|
current_point = Point(x_l93, y_l93)
|
|
df.at[i, "Longitude_Lambert"] = x_l93
|
|
df.at[i, "Latitude_Lambert"] = y_l93
|
|
|
|
previous_point = Point(wgs84_to_l93.transform(df.loc[i-1, "Longitude_WGS"],df.loc[i-1, "Latitude_WGS"]))
|
|
|
|
d1 = L.project(previous_point)
|
|
d2 = L.project(current_point)
|
|
dist = d2 - d1
|
|
df.at[i, "Distance_Cumulee"] = df.at[i-1, "Distance_Cumulee"] + dist
|
|
|
|
|
|
# === PARAMÈTRES ===
|
|
in_path = None # TODO
|
|
csv_file = "input.csv"
|
|
shapefile_name = "célé.shp"
|
|
|
|
# === 1. Chargement des données
|
|
df = pd.read_csv(os.path.join(in_path, csv_file), sep=";")
|
|
riv_line = gpd.read_file(os.path.join(in_path, shapefile_name)).to_crs(epsg=2154)
|
|
wgs84_to_l93 = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:2154", always_xy=True)
|
|
l93_to_wgs84 = pyproj.Transformer.from_crs("EPSG:2154", "EPSG:4326", always_xy=True)
|
|
|
|
# === 2. Colonnes de travail
|
|
df["has_gps"] = df["Longitude_WGS"].notna() & df["Latitude_WGS"].notna()
|
|
df["Distance_Cele"] = 0.0
|
|
|
|
# === 3. Calcul des distances entre points
|
|
project_shp(riv_line, wgs84_to_l93, l93_to_wgs84)
|
|
calcul_distance_cumulee(riv_line, wgs84_to_l93, l93_to_wgs84)
|
|
|
|
# Nettoyage
|
|
df.drop(columns=["has_gps"], inplace=True)
|
|
df.drop(columns=["Distance_entre_points_m"], inplace=True)
|
|
df.drop(columns=["Distance_Cele"], inplace=True)
|
|
|
|
|
|
# === 4. Export final
|
|
out_path = None # TODO
|
|
out_name = "output.csv"
|
|
df.to_csv(os.path.join(out_path, out_name), sep=";", index=False, encoding="utf-8")
|
|
print("✅ Fichier exporté :", out_name)
|
|
|
|
# === 5. Visualisation rapide
|
|
valid_coords = df[df["Latitude_WGS"].notna() & df["Longitude_WGS"].notna()].copy()
|
|
valid_coords["geometry"] = [Point(xy) for xy in zip(valid_coords["Longitude_WGS"], valid_coords["Latitude_WGS"])]
|
|
gdf_corr = gpd.GeoDataFrame(valid_coords, geometry="geometry", crs="EPSG:4326")
|
|
|
|
fig, ax = plt.subplots(figsize=(14, 7))
|
|
gdf_corr.plot(ax=ax, color="green", markersize=6, label="Points corrigés v2")
|
|
riv_line.to_crs(epsg=4326).plot(ax=ax, edgecolor='black', linewidth=1)
|
|
|
|
ax.set_title("Carte des points GPS corrigés projetés sur le tracé")
|
|
ax.set_xlabel("Longitude")
|
|
ax.set_ylabel("Latitude")
|
|
ax.legend()
|
|
plt.grid(True)
|
|
plt.tight_layout()
|
|
plt.show()
|