1.6) Geospatial Vector Data with geopandas

import geopandas as gpd
from shapely.geometry import Point, Polygon, LineString

stations = gpd.GeoDataFrame(
    {"name": ["Basel", "Zurich", "Bern", "Lugano"]},
    geometry=[Point(7.59, 47.56), Point(8.54, 47.37),
              Point(7.44, 46.95), Point(8.95, 46.00)],
    crs="EPSG:4326",   # WGS84 longitude/latitude in degrees
)
print(stations)
print("crs:", stations.crs.to_epsg(), "| geometry type:", stations.geom_type.iloc[0])

from pathlib import Path

Path("_files").mkdir(exist_ok=True)
stations.to_file("_files/stations.geojson", driver="GeoJSON")
stations.to_file("_files/stations.gpkg", driver="GPKG")

reread = gpd.read_file("_files/stations.gpkg")
print("reread shape:", reread.shape, "| crs preserved:", reread.crs.to_epsg())

# a shapefile round-trip, written and read the same way as GeoJSON/GeoPackage above
stations.to_file("_files/stations.shp")   # driver inferred from the .shp extension
reread_shp = gpd.read_file("_files/stations.shp")
print("reread shape:", reread_shp.shape, "| crs preserved:", reread_shp.crs.to_epsg())

# a shapefile is not one file: this write also creates _files/stations.shx, .dbf, .prj (and more)
print("shapefile sidecar files:", sorted(p.name for p in Path("_files").glob("stations.*")))

stations_lv95 = stations.to_crs(2056)   # reproject to metric Swiss coordinates

print("before:", stations.crs.to_epsg(), "->", "after:", stations_lv95.crs.to_epsg())
print("Basel WGS84 (deg):", round(stations.geometry.x.iloc[0], 3), round(stations.geometry.y.iloc[0], 3))
print("Basel LV95 (m):  ", round(stations_lv95.geometry.x.iloc[0], 1), round(stations_lv95.geometry.y.iloc[0], 1))

zones = gpd.GeoDataFrame(
    {"hazard": ["flood", "flood"]},
    geometry=[Polygon([(7, 47), (8, 47), (8, 48), (7, 48)]),      # north-west box
              Polygon([(8, 47), (9, 47), (9, 48), (8, 48)])],     # north-east box
    crs="EPSG:4326",
)

# which stations fall within a hazard zone?
in_hazard = gpd.sjoin(stations, zones, predicate="within", how="inner")
print("stations in a hazard zone:", in_hazard["name"].tolist())

# buffer: a 20 km catchment around each station (metric crs required)
buffers = stations_lv95.buffer(20_000)
print("buffer area:", round(buffers.area.iloc[0] / 1e6, 1), "km^2")

# dissolve: merge the two flood boxes into one polygon by shared attribute
dissolved = zones.dissolve(by="hazard")
print("polygons after dissolve:", len(dissolved))

# centroid: a single representative point for the dissolved polygon
zone_centroid = dissolved.geometry.iloc[0].centroid
print("dissolved hazard-zone centroid:", round(zone_centroid.x, 3), round(zone_centroid.y, 3))

# overlay: intersection of the hazard zones with the Basel buffer
basel_buffer = gpd.GeoDataFrame(geometry=[buffers.iloc[0]], crs=2056)
clip = gpd.overlay(zones.to_crs(2056), basel_buffer, how="intersection")
print("overlay produced", len(clip), "piece(s), area",
      round(clip.area.sum() / 1e6, 1), "km^2")

# length: distance measured along a line, e.g. a straight-line route between two stations
route = LineString([stations_lv95.geometry.iloc[0], stations_lv95.geometry.iloc[2]])  # Basel -> Bern
print("Basel-Bern straight-line route length:", round(route.length / 1000, 1), "km")

%matplotlib inline
import matplotlib.pyplot as plt

flagged = stations.assign(in_hazard=stations["name"].isin(in_hazard["name"]))

fig, ax = plt.subplots(figsize=(5, 5))
zones.boundary.plot(ax=ax, color="tab:red", linewidth=1)
flagged.plot(ax=ax, column="in_hazard", categorical=True, legend=True, markersize=40)
ax.set_xlabel("longitude (°E)")
ax.set_ylabel("latitude (°N)")
ax.set_title("stations and hazard zones (EPSG:4326)")
plt.show()

# the "distance between two stations" as an assistant might compute it
import warnings

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    dist_deg = stations.geometry.distance(stations.geometry.iloc[2])   # Basel..Bern etc., in EPSG:4326
    area_deg = gpd.GeoSeries(
        [Polygon([(7, 47), (8, 47), (8, 48), (7, 48)])], crs="EPSG:4326").area

print("Basel->Bern distance in degrees (wrong):", round(dist_deg.iloc[0], 4))
print("1x1 degree box area in deg^2 (wrong):", round(area_deg.iloc[0], 4))
print("geopandas warned about a geographic crs:",
      any("geographic CRS" in str(warning.message) for warning in caught))

print(f'\n{caught[0].message}')

dist_m = stations_lv95.geometry.distance(stations_lv95.geometry.iloc[2])
area_m2 = gpd.GeoSeries(
    [Polygon([(7, 47), (8, 47), (8, 48), (7, 48)])], crs="EPSG:4326").to_crs(2056).area

print("Basel->Bern distance:", round(dist_m.iloc[0] / 1000, 1), "km")
print("1x1 degree box area:", round(area_m2.iloc[0] / 1e6, 0), "km^2")