Código
import pandas as pd
import geopandas as gpd
import folium
from shapely.geometry import Polygon
from pathlib import Path
# Input CSV with geographic coordinates (Longitude/Latitude in EPSG:4326)
csv_path = "field_data/rugu_arbeyal.csv"
cover_csv = "field_data/rugu_arbeyal_cover.csv"
# Load tabular data
df = pd.read_csv(csv_path)
df_cover = pd.read_csv(cover_csv)
df_cover["cover_pct"] = pd.to_numeric(df_cover["Description"], errors="coerce")
df_cover = df_cover.dropna(subset=["cover_pct", "Longitude", "Latitude"])
df_cover = df_cover[df_cover["cover_pct"].between(0, 100, inclusive="both")].copy()
# Build point GeoDataFrames
gdf_pts = gpd.GeoDataFrame(
df.copy(),
geometry=gpd.points_from_xy(df["Longitude"], df["Latitude"]),
crs="EPSG:4326",
)
gdf_cover = gpd.GeoDataFrame(
df_cover.copy(),
geometry=gpd.points_from_xy(df_cover["Longitude"], df_cover["Latitude"]),
crs="EPSG:4326",
)
# Build two polygons from point IDs
gdf_pts["point_id"] = pd.to_numeric(gdf_pts["Name"], errors="coerce")
def build_polygon(point_min, point_max):
part = gdf_pts[
(gdf_pts["point_id"] >= point_min) & (gdf_pts["point_id"] <= point_max)
].sort_values("point_id")
coords = list(zip(part["Longitude"], part["Latitude"]))
if coords and coords[0] != coords[-1]:
coords.append(coords[0])
return Polygon(coords)
poly_1 = build_polygon(1, 27)
poly_2 = build_polygon(28, 127)
gdf_poly = gpd.GeoDataFrame(
{"name": ["polygon_1", "polygon_2"]},
geometry=[poly_1, poly_2],
crs="EPSG:4326",
)
# Compute total area in square meters using projected CRS
total_area_m2 = gdf_poly.to_crs(epsg=25830).area.sum()
# Save polygons to GeoJSON in results folder
results_dir = Path("results")
results_dir.mkdir(parents=True, exist_ok=True)
geojson_out = results_dir / "rugu_arbeyal_polygons.geojson"
gdf_poly.to_file(geojson_out, driver="GeoJSON")
# Build interactive Folium map centered on polygon centroid
center = gdf_poly.geometry.iloc[0].centroid
m = folium.Map(location=[center.y, center.x], zoom_start=18, tiles=None)
folium.TileLayer(
tiles="https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
attr="Tiles © Esri",
name="Esri World Imagery",
overlay=False,
control=True,
).add_to(m)
# Add polygon
folium.GeoJson(
gdf_poly.__geo_interface__,
name="Boundary polygon",
style_function=lambda _feature: {
"color": "red",
"weight": 3,
"fillColor": "red",
"fillOpacity": 0.1,
},
).add_to(m)
# Add boundary points
for row in gdf_pts.itertuples():
point_name = getattr(row, "Name", "Unknown")
folium.CircleMarker(
location=[row.geometry.y, row.geometry.x],
radius=4,
color="black",
weight=1,
fill=True,
fill_color="yellow",
fill_opacity=1.0,
popup=f"Name: {point_name}",
tooltip=str(point_name),
).add_to(m)
# Add % cover points
for row in gdf_cover.itertuples():
folium.CircleMarker(
location=[row.geometry.y, row.geometry.x],
radius=4,
color="black",
weight=1,
fill=True,
fill_color="#2ca25f",
fill_opacity=1.0,
popup=f"Cobertura: {row.cover_pct:.1f}%",
tooltip=f"{row.cover_pct:.1f}%",
).add_to(m)
# Zoom to polygon extent so all vertices are visible
minx, miny, maxx, maxy = gdf_poly.total_bounds
m.fit_bounds([[miny, minx], [maxy, maxx]])
# Legend with red line symbol and total area label
legend_html = f"""
<div style="
position: fixed;
bottom: 20px;
left: 20px;
z-index: 1000;
background-color: white;
border: 1px solid #666;
border-radius: 6px;
padding: 10px 12px;
font-size: 14px;
box-shadow: 0 1px 6px rgba(0, 0, 0, 0.25);
">
<div style="display: flex; align-items: center; gap: 8px;">
<span style="display: inline-block; width: 26px; border-top: 3px solid red;"></span>
<span><b>Área arribazón <br>{total_area_m2:,.2f} m2</b></span>
</div>
<div style="display: flex; align-items: center; gap: 8px; margin-top: 8px;">
<span style="display: inline-block; width: 9px; height: 9px; border-radius: 50%; background: yellow; border: 1px solid black;"></span>
<span><b>Bordes</b></span>
</div>
<div style="display: flex; align-items: center; gap: 8px; margin-top: 8px;">
<span style="display: inline-block; width: 9px; height: 9px; border-radius: 50%; background: #2ca25f; border: 1px solid black;"></span>
<span><b>% cobertura</b></span>
</div>
</div>
"""
m.get_root().html.add_child(folium.Element(legend_html))
# Show map in Quarto
mMake this Notebook Trusted to load map: File -> Trust Notebook