python 地理信息 python 地形

写在前面

Python使用Cartopy绘制全球地形图

效果图

 

1 不含中国边界和省界的全球地形图

直接参考公众号的推送就没问题

2 含中国边界和省界的中国地形图

（综合1和之前画中国站点的代码就可以）

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib import rcParams
from cartopy.io import shapereader as shpreader


config = {"font.family": 'Times New Roman',
          "font.size": 14, "mathtext.fontset": 'stix'}
rcParams.update(config)

# 读取全球地形数据
ds = xr.open_dataset("D:桌面\ETOPO2v2c_f4.nc")

# 准备用于绘图的数据
lon = np.linspace(min(ds['x'].data), max(
    ds['x'].data), len(ds['x'].data))  # 经度
lat = np.linspace(min(ds['y'].data), max(
    ds['y'].data), len(ds['y'].data))  # 纬度
dem = ds['z'].data

# 构建经纬网格
lon, lat = np.meshgrid(lon, lat)

# 设置地图全局属性
fig = plt.figure(1, figsize=[16, 9])
proj = ccrs.PlateCarree()
ax = plt.subplot(1, 1, 1, projection=proj)
extent = [70, 136, 15, 60]
ax.set_extent(extent, crs=proj)

china = shpreader.Reader(r"D:bou2_4l.dbf").geometries()
# 绘制中国国界省界九段线等等


# 创建底图，设置地图投影为World Plate Carrée，分辨率为高分辨率，地图范围为全球 20个
levels = [-8000, -6000, -4000, -2000, -1000, -200, -50, 0, 50,
          200, 500, 1000, 1500, 2000, 3000, 4000, 5000, 6000, 7000, 8000]
# levels=np.arange(-8000,8000,1000)

# 创建分级
color = ['#084594', '#2171b5', '#4292c6', '#6baed6', '#9ecae1', '#c6dbef', '#deebf7', '#006837',
         '#31a354', '#78c679', '#addd8e', '#d9f0a3', '#f7fcb9', '#c9bc87', '#a69165', '#856b49',
         '#664830', '#ad9591', '#d7ccca']  # 设置色带，19个颜色
cf = ax.contourf(lon, lat, dem, levels=levels, colors=color, extend='both')


ax.add_geometries(china, ccrs.PlateCarree(),
                  facecolor='none', edgecolor='black', zorder=1)  # , linewidth=12)  # , zorder=1)

# 设置图例,shrink调整色标长度
cb = plt.colorbar(cf, shrink=0.7, orientation='horizontal',
                  pad=0.05, aspect=50, ticks=levels, extend='both')
cb.set_ticks(levels)  # 设置色带刻度
cb.ax.tick_params(labelsize=10)
cb.set_label('Global Elevation(meter)', fontsize=12)  # 设置图例名称和字体大小
ax.set_xticks(np.arange(extent[0], extent[1]+1, 10), crs=proj)
ax.set_yticks(np.arange(extent[-2], extent[-1]+1, 10), crs=proj)
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
ax.yaxis.set_major_formatter(LatitudeFormatter())
plt.savefig('D:\何卓尧\桌面\plot158.png', dpi=300,
            bbox_inches='tight', pad_inches=0)
plt.show()

注意一定要先contourf再画边界，不然边界显现不出来。

add_geometries里facecolor要设为none，可以尝试一下设为black，会导致边界不是均匀的线。