如何用Python沿海岸线每隔100英里绘制里程桩？

要实现沿海岸线每隔100英里绘制里程桩，同时给桩点间海岸线着色的需求，核心是将海岸线的地理线段转换为带累计距离的路径，按固定间隔拆分路径后分别绘制分段线条和里程桩。具体步骤如下：

1. 获取并预处理海岸线数据 #

Cartopy自带的海岸线数据来自Natural Earth，可通过cartopy.feature.NaturalEarthFeature直接调用，或用cartopy.io.shapereader读取shp文件。获取后将数据转换为Shapely的LineString/MultiLineString对象，方便后续距离计算。

2. 按固定间隔拆分海岸线路径 #

地理坐标下的距离需用大圆距离计算实际地表长度，避免欧氏距离的误差。通过累加相邻点位的实际距离，筛选出每隔100英里的拆分点，同时将原始海岸线拆分为对应间隔的子线段。

3. 分段着色海岸线并绘制里程桩 #

给拆分后的子线段赋予不同颜色（可自定义着色规则），再在拆分点位置绘制标记和里程数标注。

完整代码示例 #

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import numpy as np
from shapely.geometry import LineString, MultiLineString
from geopy.distance import geodesic

# 单位转换：1英里 = 1609.34米
MILE_TO_METER = 1609.34
INTERVAL_MILE = 100
INTERVAL_METER = INTERVAL_MILE * MILE_TO_METER

def split_line_by_distance(line, interval):
    """将LineString按固定距离间隔拆分，返回拆分点坐标和分段线段"""
    coords = list(line.coords)
    if len(coords) < 2:
        return [], []
    
    split_points = [coords[0]]
    segments = []
    current_dist = 0
    prev_coord = coords[0]
    
    for coord in coords[1:]:
        # 计算当前点与上一点的大圆距离（米）
        dist = geodesic(prev_coord, coord).meters
        if current_dist + dist >= interval:
            # 计算拆分位置的坐标
            ratio = (interval - current_dist) / dist
            split_lon = prev_coord[0] + ratio * (coord[0] - prev_coord[0])
            split_lat = prev_coord[1] + ratio * (coord[1] - prev_coord[1])
            split_point = (split_lon, split_lat)
            split_points.append(split_point)
            # 添加分段线段
            segments.append(LineString([prev_coord, split_point]))
            # 更新累计距离和基准点
            current_dist = current_dist + dist - interval
            prev_coord = split_point
        else:
            current_dist += dist
            segments.append(LineString([prev_coord, coord]))
            prev_coord = coord
    
    # 补充最后一个端点
    if prev_coord != coords[-1]:
        split_points.append(coords[-1])
    return split_points, segments

# 设置地图投影（示例为经纬度投影，可按需调整）
proj = ccrs.PlateCarree()
fig, ax = plt.subplots(figsize=(12, 8), subplot_kw={'projection': proj})

# 添加基础海岸线（浅色细线条）
coastline_base = cfeature.NaturalEarthFeature(
    category='physical',
    name='coastline',
    scale='10m',
    edgecolor='#cccccc',
    facecolor='none'
)
ax.add_feature(coastline_base, linewidth=0.5)

# 读取海岸线几何数据
reader = Reader(cfeature.NaturalEarthFeature('physical', 'coastline', '10m').path)
geoms = list(reader.geometries())

# 准备渐变颜色列表
colors = plt.cm.viridis(np.linspace(0, 1, 20))
color_idx = 0
all_mile_markers = []

# 处理每一段海岸线几何对象
for geom in geoms:
    lines = list(geom.geoms) if isinstance(geom, MultiLineString) else [geom]
    for line in lines:
        split_points, segments = split_line_by_distance(line, INTERVAL_METER)
        if len(split_points) < 2:
            continue
        
        # 绘制分段着色的海岸线
        for seg in segments:
            ax.plot(*seg.xy, transform=proj, color=colors[color_idx % len(colors)], linewidth=2)
            color_idx += 1
        
        # 收集里程桩点位（排除首尾重复点）
        all_mile_markers.extend(split_points[1:-1])

# 绘制里程桩标记与文字标注
for i, (lon, lat) in enumerate(all_mile_markers):
    ax.scatter(lon, lat, transform=proj, color='red', s=50, zorder=5)
    ax.text(lon + 0.15, lat + 0.15, f"{(i+1)*INTERVAL_MILE}mi", 
            transform=proj, fontsize=8, zorder=6, bbox={'facecolor':'white', 'alpha':0.7, 'pad':2})

# 设置地图显示范围（示例为北美东海岸，可自定义）
ax.set_extent([-80, -60, 30, 45], crs=proj)
ax.gridlines(draw_labels=True, linestyle='--', alpha=0.7)

plt.show()

关键说明 #

• split_line_by_distance是核心函数：通过geodesic计算实际地表距离，确保间隔的准确性，避免地理坐标下的欧氏距离误差。
• 投影选择：示例用经纬度投影，若处理大范围区域，建议改用等距投影（如UTM）提升计算效率，需根据目标区域调整投影参数。
• 着色逻辑：示例用渐变颜色循环，可根据需求修改为固定颜色、按里程区间着色等规则。

内容的提问来源于stack exchange，提问作者Feng Hu