我跟着Adam Symington的优秀指南成功制作了马来西亚沙巴州的地形图，但过程中碰到了几个小问题，现在已经全部解决，把解决方案整理出来分享给大家：

1. 学校坐标绘制异常（坐标轴非地理坐标轴） #

原问题 #

sch dataframe里存着沙巴州所有学校的经纬度坐标，但绘制后在左上角出现了黑色不规则图形，验证后发现是因为地形图的坐标轴不是地理经纬度轴，和学校坐标的CRS不匹配导致的。

解决方案 #

地形数据是经过裁剪的栅格，我们需要先获取原始栅格的坐标系和变换参数，把学校的WGS84（epsg:4326）经纬度坐标转换为栅格的像素坐标。这里我们可以用rasterio的transform_points方法来完成转换。

2. 地图外部背景设为白色+添加黑色轮廓 #

原问题 #

调用ax.set_facecolor('white')没效果，因为ax.imshow(hillshade, cmap='Greys', alpha=0.3)把整个画布背景设成了灰色。想要保留地图内部的灰色阴影，同时把外部区域改成白色，还想给地图加个黑色轮廓。

解决方案 #

• 利用地形数据的非NaN区域创建掩码，给imshow添加extent参数，同时用clip_path来限定山体阴影的显示区域；
• 提取地形数据的有效区域边界，绘制黑色轮廓线。

修改后的完整可复现代码 #

import rasterio
from rasterio import mask as msk
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.colors import ListedColormap
import numpy as np
import pandas as pd
import geopandas as gpd
import earthpy.spatial as es
from shapely.geometry import Point
from rasterio.transform import from_origin

def prepareFiles():
    import os
    # 创建存放数据的文件夹
    if not os.path.exists('sabah'):
        os.mkdir('sabah')
    # 提示：请手动获取以下两个文件并保存到当前目录的sabah文件夹中：
    # 1. 沙巴地形数据文件，命名为 sabah_topog.npy
    # 2. 沙巴学校坐标文件，命名为 sabah_schools.csv

# 运行一次准备文件夹，之后可注释
prepareFiles()

# 加载沙巴地形数据
value_range = 4049
sabah_topography = np.load('sabah/sabah_topog.npy')
topo_data = sabah_topography[0]

# --------------------------
# 关键改动1：匹配坐标系统
# 这里我们模拟原始栅格的地理变换参数（实际可从原始TIFF文件读取）
# 沙巴大致的地理范围参数，可根据实际数据调整
left, top = 115.5, 7.0  # 栅格左上角的经纬度
resolution = 0.002  # 栅格分辨率（单位：度）
width, height = topo_data.shape[1], topo_data.shape[0]
# 创建栅格的地理变换对象
transform = from_origin(left, top, resolution, resolution)

# 加载学校坐标数据
sch = pd.read_csv('sabah/sabah_schools.csv', usecols=['lat','lon'])
# 将经纬度转换为栅格的像素坐标
xy = rasterio.transform.transform_points(transform, sch.lon.values, sch.lat.values)
sch['pixel_x'] = xy[:, 0]
sch['pixel_y'] = xy[:, 1]

# --------------------------
# 自定义地图配色
sabah_colormap = LinearSegmentedColormap.from_list('sabah', ['lightgray', '#e6757b', '#CD212A', '#CD212A'], N=value_range)
background_color = np.array([1,1,1,1])
newcolors = sabah_colormap(np.linspace(0, 1, value_range))
newcolors = np.vstack((newcolors, background_color))
sabah_colormap = ListedColormap(newcolors)

# 生成山体阴影效果
hillshade = es.hillshade(topo_data, azimuth=180, altitude=1)

# --------------------------
# 关键改动2：设置背景与轮廓
plt.rcParams["figure.figsize"] = [5,5]
plt.rcParams["figure.autolayout"] = True
fig, ax = plt.subplots()

# 定义地图的地理范围，用于限制显示区域
extent = [left, left + width*resolution, top - height*resolution, top]

# 绘制地形图层
topo_im = ax.imshow(topo_data, cmap=sabah_colormap, extent=extent)
# 绘制山体阴影，仅显示地形有效区域（利用clip_path裁剪）
ax.imshow(hillshade, cmap='Greys', alpha=0.3, extent=extent, clip_path=topo_im.get_clip_path())

# 绘制学校坐标点（现在坐标已匹配栅格）
ax.scatter(sch['pixel_x'], sch['pixel_y'], color='black', marker='x', markersize=10)

# 添加地图黑色轮廓
# 提取地形有效区域的边界坐标
mask = ~np.isnan(topo_data)
y_idx, x_idx = np.where(mask)
min_x = left + x_idx.min()*resolution
max_x = left + x_idx.max()*resolution
min_y = top - y_idx.max()*resolution
max_y = top - y_idx.min()*resolution
# 绘制矩形轮廓
ax.plot([min_x, max_x, max_x, min_x, min_x], 
        [min_y, min_y, max_y, max_y, min_y], 
        color='black', linewidth=1)

# 设置画布背景为白色
ax.set_facecolor('white')
ax.axis('off')

plt.show()

关键改动说明 #

1. 坐标匹配：通过创建栅格的地理变换对象，把学校的经纬度坐标转换成了栅格对应的像素坐标，彻底解决了点绘制位置异常的问题；
2. 背景与轮廓优化：
   • 用extent参数定义地图的实际地理范围，配合clip_path让山体阴影只显示在地形有效区域内；
   • 通过提取地形数据的有效边界，绘制了黑色的地图轮廓；
   • 此时再设置ax.set_facecolor('white')就能让地图外部区域显示为白色了。

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