基于QGIS与QuickMapServices模拟海平面上升生成海岸线轮廓

可以实现，但需注意：OpenTopoMap是渲染后的栅格瓦片（图片格式），而非原始DEM（数字高程模型）数据，因此模拟海平面上升的精度会受图像分析的限制，仅能生成近似结果。如果需要高精度的海平面上升模拟，建议配合DEM数据，但如果坚持仅用OpenTopoMap，以下是Python程序化实现步骤：

1. 准备依赖库 #

先安装所需工具包：

pip install requests rasterio geopandas matplotlib shapely pillow opencv-python svgwrite pyproj

2. 获取目标区域边界（WA、OR、BC） #

用geopandas加载行政区划数据，裁剪出美国华盛顿州、俄勒冈州及加拿大不列颠哥伦比亚省的范围：

import geopandas as gpd
from pyproj import Transformer

# 加载全球低分辨率行政区划数据
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

# 筛选目标区域
wa_or_bc = world[
    ((world['admin'] == 'United States of America') & (world['subunit'].isin(['Washington', 'Oregon']))) |
    ((world['admin'] == 'Canada') & (world['subunit'] == 'British Columbia'))
]

# 转换为OpenTopoMap使用的Web Mercator投影（EPSG:3857）
wa_or_bc = wa_or_bc.to_crs(epsg=3857)
# 获取区域边界框
bbox = wa_or_bc.total_bounds

3. 下载并拼接OpenTopoMap瓦片 #

根据目标区域的边界计算所需瓦片范围，下载后拼接成完整栅格图像：

import requests
from PIL import Image
from io import BytesIO
import math
import numpy as np

def deg2num(lat_deg, lon_deg, zoom):
    """经纬度转瓦片坐标"""
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
    return (xtile, ytile)

# 将Web Mercator边界转为经纬度，用于计算瓦片
transformer = Transformer.from_crs(3857, 4326, always_xy=True)
min_lon, min_lat = transformer.transform(bbox[0], bbox[1])
max_lon, max_lat = transformer.transform(bbox[2], bbox[3])

# 选择瓦片缩放级别（11级平衡精度与文件大小，可调整）
zoom = 11
min_x, max_y = deg2num(max_lat, min_lon, zoom)
max_x, min_y = deg2num(min_lat, max_lon, zoom)

# 下载所有瓦片
tiles = []
for x in range(min_x, max_x + 1):
    row = []
    for y in range(min_y, max_y + 1):
        url = f"https://tile.opentopomap.org/{zoom}/{x}/{y}.png"
        response = requests.get(url)
        img = Image.open(BytesIO(response.content))
        row.append(img)
    tiles.append(row)

# 拼接瓦片为完整图像
tile_w, tile_h = tiles[0][0].size
full_w, full_h = tile_w * len(tiles[0]), tile_h * len(tiles)
full_img = Image.new('RGB', (full_w, full_h))

for i, row in enumerate(tiles):
    for j, tile in enumerate(row):
        full_img.paste(tile, (j * tile_w, i * tile_h))

# 裁剪到目标区域边界
from rasterio.transform import from_bounds
from rasterio.features import geometry_mask

transform = from_bounds(bbox[0], bbox[1], bbox[2], bbox[3], full_w, full_h)
clip_poly = wa_or_bc.unary_union
mask = geometry_mask([clip_poly], out_shape=(full_h, full_w), transform=transform, invert=True)

full_img_np = np.array(full_img)
full_img_np[~mask] = 255  # 非目标区域设为白色
clipped_img = Image.fromarray(full_img_np)

4. 模拟海平面上升并提取海岸线 #

通过图像边缘检测提取海陆交界线，再根据海平面上升高度计算像素偏移量，生成新海岸线：

import cv2

# 转换为灰度图并提取边缘
gray = cv2.cvtColor(np.array(clipped_img), cv2.COLOR_RGB2GRAY)
edges = cv2.Canny(gray, 50, 150)
# 去除噪声
kernel = np.ones((3,3), np.uint8)
edges = cv2.morphologyEx(edges, cv2.MORPH_CLOSE, kernel)

# 计算米/像素转换比例（Web Mercator投影下）
lat_center = (min_lat + max_lat) / 2
meters_per_pixel = (40075016.686 * math.cos(math.radians(lat_center))) / (2**zoom * 256)
x = 10  # 替换为你需要模拟的海平面上升米数
offset_pixels = int(x / meters_per_pixel)

# 生成海平面上升后的海岸线（向内偏移）
offset_kernel = np.ones((offset_pixels*2, offset_pixels*2), np.uint8)
flooded_area = cv2.erode(255 - edges, offset_kernel)
new_coastline = cv2.Canny(flooded_area, 50, 150)

5. 生成4K分辨率黑色线条图 #

将提取的海岸线转换为4K分辨率的PNG或SVG格式：

# 生成PNG
result_img = np.full((new_coastline.shape[0], new_coastline.shape[1], 3), 255, dtype=np.uint8)
result_img[new_coastline > 0] = [0, 0, 0]  # 线条设为黑色
# 调整为4K分辨率（3840x2160）
result_img = cv2.resize(result_img, (3840, 2160), interpolation=cv2.INTER_LINEAR)
cv2.imwrite('coastline_4k.png', result_img)

# 生成SVG矢量图
import svgwrite

dwg = svgwrite.Drawing('coastline_4k.svg', size=(3840, 2160))
contours, _ = cv2.findContours(new_coastline, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
for contour in contours:
    path_data = "M " + " L ".join([f"{p[0][0]} {p[0][1]}" for p in contour]) + " Z"
    dwg.add(dwg.path(d=path_data, fill='none', stroke='black', stroke_width=1))
dwg.save()

关键提示 #

• OpenTopoMap的渲染颜色可能因区域、缩放级别略有差异，需根据实际效果调整边缘检测的阈值参数。
• 若需要更高精度的海平面模拟，建议配合DEM数据（如NASA SRTM）提取等高线，再叠加到OpenTopoMap上。

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