如何基于Basemap对经纬度散点数据应用高斯滤波绘制平滑地图

你要实现论文中的300km高斯滤波、1km分辨率平滑地图效果，不要直接套用普通位图的高斯滤波方案，核心流程为「投影坐标转换→目标分辨率规则网格插值→匹配实际地理距离的高斯卷积→栅格图层绘制」，具体实现如下：

提前安装需要的库：numpy、scipy、matplotlib、mpl_toolkits.basemap

1. 初始化Basemap并转换投影坐标 #

经纬度是球面坐标，无法直接按平面距离计算滤波参数，先把所有散点转换为你所用LCC投影下的平面坐标（单位：米），保证距离计算准确：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
from scipy.ndimage import gaussian_filter
from mpl_toolkits.basemap import Basemap

# 原有研究区范围参数
la=35.12
ua=35.6
ln=127.1
un=128
m_lat=(la+ua)/2
m_lon=(ln+un)/2

# 初始化Basemap
m = Basemap(
    projection='lcc', 
    lat_0=m_lat, 
    lon_0=m_lon, 
    resolution='f', 
    llcrnrlon=ln, 
    llcrnrlat=la, 
    urcrnrlon=un, 
    urcrnrlat=ua
)

# 把散点经纬度转成投影平面坐标
x_proj, y_proj = m(X_list, Y_list)

2. 生成1km分辨率的规则网格 #

按需要的1km空间分辨率，生成覆盖整个研究区的规则网格，网格间隔对应1000米：

# 获取研究区在投影坐标系下的边界
x_min, y_min = m(ln, la)
x_max, y_max = m(un, ua)

# 生成1km间隔的网格矩阵
grid_x, grid_y = np.meshgrid(
    np.arange(x_min, x_max, 1000),  # 1000米=1km，修改该值可调整输出分辨率
    np.arange(y_min, y_max, 1000)
)

3. 散点线性插值到网格 #

对应论文提到的线性插值步骤，把离散点的Z值插值到规则网格上，边缘空值用最近邻插值补全，避免出现空洞：

# 先做线性插值
grid_z = griddata(
    points=(x_proj, y_proj),
    values=Z_list,
    xi=(grid_x, grid_y),
    method='linear'
)

# 用最近邻插值补线性插值产生的边缘空值
grid_z_nearest = griddata(
    points=(x_proj, y_proj),
    values=Z_list,
    xi=(grid_x, grid_y),
    method='nearest'
)
grid_z[np.isnan(grid_z)] = grid_z_nearest[np.isnan(grid_z)]

4. 300km全宽高斯滤波 #

注意论文提到的300km是高斯核的全宽半高（FWHM），需要先换算成高斯滤波的sigma参数，再结合网格分辨率换算为像素单位的sigma值做卷积：

# 参数配置
filter_fwhm_km = 300  # 高斯滤波全宽，和论文参数对齐
grid_res_km = 1       # 网格空间分辨率
# FWHM转sigma的固定换算公式：sigma = FWHM / (2*sqrt(2*ln2)) ≈ FWHM/2.355
sigma_pixel = (filter_fwhm_km / 2.355) / grid_res_km

# 做高斯平滑，truncate=4表示高斯核截断在4倍sigma位置，足够覆盖有效滤波范围
grid_z_smooth = gaussian_filter(grid_z, sigma=sigma_pixel, truncate=4, mode='nearest')

5. 绘制平滑地图 #

替换原来的scatter散点绘制代码，用滤波后的网格数据绘制平滑面图层：

# 按需绘制底图要素
m.drawcoastlines(linewidth=0.8)
m.drawcountries(linewidth=0.8)

# 绘制平滑后的属性面
im = m.pcolormesh(
    grid_x, grid_y, grid_z_smooth,
    cmap='bwr',
    vmin=-1,
    vmax=1,
    shading='auto'
)

# 添加色带
plt.colorbar(im, label='属性值Z')
plt.show()

• 调整输出地图分辨率：修改生成网格时的步长值即可，比如要500m分辨率就把步长从1000改成500，同步重新计算sigma_pixel即可
• 调整滤波平滑程度：修改filter_fwhm_km数值，数值越大平滑效果越强
• 如果边缘出现异常值，可以在gaussian_filter里调整mode参数，可选reflect、wrap等模式适配数据分布

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