Python地图投影至lon=30时格陵兰出现异常线条的技术问询

问题原因 #

格陵兰上空的异常线条是跨投影断点的多边形在转换时产生的连接错误。当投影中心设为30°E时，投影的“边界”落在30°±180°（即-150°W和150°E），格陵兰的几何形状跨越了这个边界，导致GeoPandas在投影转换时生成错误的连线。

解决方案 #

以下是两种可行的修复方法，结合你的代码修改：

方法1：几何修复+投影后清理

在投影前后对几何数据进行有效性修复，过滤空几何：

import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import seaborn as sns
import geopandas as gpd
import os
import zipfile
import pyproj

# 定义A1纸张尺寸（英寸）
a1_width_inches = 33.1
a1_height_inches = 23.4

# 指定shapefile的URL和文件名
shapefile_url = "https://naturalearth.s3.amazonaws.com/110m_cultural/ne_110m_admin_0_countries.zip"
shapefile_zip = "ne_110m_admin_0_countries.zip"
shapefile_name = "ne_110m_admin_0_countries.shp"

# 如果shapefile不存在则下载
if not os.path.exists(shapefile_name):
    # 下载压缩文件
    !wget {shapefile_url} -O {shapefile_zip}

    # 解压shapefile
    with zipfile.ZipFile(shapefile_zip, 'r') as zip_ref:
        zip_ref.extractall()

    print(f"Shapefile已下载并解压至: {shapefile_name}")

# 从下载的文件中加载世界地图数据
world = gpd.read_file(shapefile_name)

# 修复原始几何的有效性问题，过滤空几何
world = world.make_valid()
world = world[~world.is_empty]

# 定义温克尔三重投影（Winkel Tripel）
winkel_tripel_crs = "+proj=wintri +lon_0=30 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

# 对GeoDataFrame进行投影转换
world = world.to_crs(winkel_tripel_crs)

# 投影后再次修复几何，去除异常连接
world = world.make_valid()
world = world[~world.is_empty]

方法2：平移经度避免跨投影断点

在投影前将所有经度平移，让投影中心（30°E）成为坐标中心，避免多边形跨越投影边界：

import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import seaborn as sns
import geopandas as gpd
import os
import zipfile
import pyproj
from shapely.ops import transform

# 定义A1纸张尺寸（英寸）
a1_width_inches = 33.1
a1_height_inches = 23.4

# 指定shapefile的URL和文件名
shapefile_url = "https://naturalearth.s3.amazonaws.com/110m_cultural/ne_110m_admin_0_countries.zip"
shapefile_zip = "ne_110m_admin_0_countries.zip"
shapefile_name = "ne_110m_admin_0_countries.shp"

# 如果shapefile不存在则下载
if not os.path.exists(shapefile_name):
    # 下载压缩文件
    !wget {shapefile_url} -O {shapefile_zip}

    # 解压shapefile
    with zipfile.ZipFile(shapefile_zip, 'r') as zip_ref:
        zip_ref.extractall()

    print(f"Shapefile已下载并解压至: {shapefile_name}")

# 从下载的文件中加载世界地图数据
world = gpd.read_file(shapefile_name)

# 定义经度平移函数：将坐标中心移至30°E，避免跨投影断点
def shift_to_center(geom, center_lon=30):
    def shift_coords(x, y, z=None):
        x_shifted = x + center_lon
        if x_shifted > 180:
            x_shifted -= 360
        elif x_shifted < -180:
            x_shifted += 360
        return (x_shifted, y)
    return transform(shift_coords, geom)

# 应用平移，重新设置CRS
world['geometry'] = world['geometry'].apply(shift_to_center)
world = world.set_crs("EPSG:4326")

# 定义温克尔三重投影（Winkel Tripel）
winkel_tripel_crs = "+proj=wintri +lon_0=30 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

# 对GeoDataFrame进行投影转换
world = world.to_crs(winkel_tripel_crs)

说明 #

• 方法1适合快速修复，通过make_valid()处理几何错误，过滤空几何来消除异常线条；
• 方法2从根源避免跨投影边界的问题，平移经度后再投影，能更彻底解决此类问题。

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