绘制带小地球的海洋地形图

在本篇博客中，我们将介绍如何使用 Python 绘制一个包含全球海洋地形的地图，并通过 Cartopy 和 GDAL 库展示海底水深数据。为了更好地理解海洋地理信息，我们还将嵌入一个小地球的 「Orthographic 投影」，以突出显示感兴趣的区域。提高科研绘图的高级感。

安装所需库

首先，确保你已经安装了以下依赖：

pip install numpy matplotlib cartopy gdal shapely

步骤 1: 导入必要的库

我们将使用 NumPy 进行数据处理，Matplotlib 用于绘图，Cartopy 用于地图投影和海洋特征绘制，GDAL 用于读取和处理地形数据。

import numpy as npimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom osgeo import osr, gdalfrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterfrom mpl_toolkits.axes_grid1.inset_locator import inset_axesfrom shapely.geometry import boxfrom cartopy.feature import ShapelyFeaturefrom cartopy.mpl.geoaxes import GeoAxes

步骤 2: 设置地图显示范围和随机散点数据

我们定义了地图的显示范围，并生成了一些随机的散点数据来演示如何在地图上标记特定位置。

# 设置绘图范围lonmin = 106lonmax = 124latmin = 12latmax = 25# 随机散点数据（演示用）np.random.seed(42)num_points = 100longitudes = np.random.uniform(lonmin + 5, lonmax - 4, num_points)latitudes = np.random.uniform(latmin + 4, latmax - 4, num_points)angles = np.random.uniform(0, 9, num_points)

步骤 3: 读取并处理 ETOPO 地形数据

接下来，我们读取全球海底地形数据集 「ETOPO1」，并裁剪出指定区域的水深数据。

# 读取 ETOPO 地形数据gdal.AllRegister()orig_filepath = r'ETOPO1_Bed_g_geotiff.tif'ds = gdal.Open(orig_filepath)band = ds.GetRasterBand(1)imgarr = band.ReadAsArray()width = ds.RasterXSizeheight = ds.RasterYSizegt = ds.GetGeoTransform()proj = ds.GetProjection()lon_min = gt[0]lat_min = gt[3] + width * gt[4] + height * gt[5]lon_max = gt[0] + width * gt[1] + height * gt[2]lat_max = gt[3]pcs = osr.SpatialReference()pcs.ImportFromWkt(proj)lon_ori = np.linspace(lon_min, lon_max, width)lat_ori = np.linspace(lat_min, lat_max, height)lons_ori, lats_ori = np.meshgrid(lon_ori, lat_ori[::-1])imgarr = np.array(imgarr)

步骤 4: 剪裁并处理水深数据

我们根据设定的经纬度范围裁剪出目标区域的水深数据，并将陆地区域的值设置为零。

# 裁剪水深latloc = [(90 - latmax) * 60, (90 - latmin) * 60]lonloc = [(180 + lonmin) * 60, (180 + lonmax) * 60]waterdepth = imgarr[int(latloc[0]):int(latloc[1]), int(lonloc[0]):int(lonloc[1])]waterdepth[waterdepth >= 0] = 0lon1 = lons_ori[int(latloc[0]):int(latloc[1]), int(lonloc[0]):int(lonloc[1])]lat1 = lats_ori[int(latloc[0]):int(latloc[1]), int(lonloc[0]):int(lonloc[1])]land_mask = waterdepth >= 0

步骤 5: 绘制主地图和小地球

现在，我们绘制主地图，添加海岸线、陆地和网格线，并绘制小地球（Orthographic 投影），显示感兴趣区域。

# 主图fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={'projection': ccrs.PlateCarree()})plt.rcParams['font.sans-serif'] = ['Arial']plt.rcParams['font.size'] = 18ax.set_extent([lonmin, lonmax, latmin, latmax], crs=ccrs.PlateCarree())ax.add_feature(cfeature.COASTLINE, linewidth=1.0)ax.add_feature(cfeature.LAND, color='lightgray')# 增加坐标轴刻度ax.set_xticks(np.arange(lonmin, lonmax + 1, 2), crs=ccrs.PlateCarree())ax.set_yticks(np.arange(latmin, latmax + 1, 2), crs=ccrs.PlateCarree())lon_formatter = LongitudeFormatter(zero_direction_label=True)lat_formatter = LatitudeFormatter()ax.xaxis.set_major_formatter(lon_formatter)ax.yaxis.set_major_formatter(lat_formatter)ax.tick_params(axis='both', which='major', labelsize=14,                direction='out',                    length=3,                           width=4.5)           gl = ax.gridlines(draw_labels=False, linewidth=0.1, color='gray', alpha=0.5, linestyle='--')# 海底地形图层pcm = ax.pcolormesh(lon1, lat1, waterdepth, cmap='terrain', shading='auto')pcm.set_array(np.ma.masked_array(waterdepth, mask=land_mask))

步骤 6: 添加小地球并突出显示区域

在主图上方，我们嵌入一个小地球，并通过红框标出当前展示区域。

# 小地球 insetax_globe = inset_axes(    ax,    width="150%", height="150%",    bbox_to_anchor=[-0.04, 0.82, 0.25, 0.25],    bbox_transform=ax.transAxes,    axes_class=GeoAxes,    axes_kwargs={'projection': ccrs.Orthographic(120, 20)})ax_globe.stock_img()ax_globe.set_global()# 隐藏小地球的边框和刻度for spine in ax_globe.spines.values():    spine.set_visible(False)ax_globe.set_xticks([])ax_globe.set_yticks([])# 在小地球上加红框，标出主图区域bbox = box(lonmin, latmin, lonmax, latmax)feature = ShapelyFeature([bbox], ccrs.PlateCarree(), edgecolor='red', facecolor='none', linewidth=1.2)ax_globe.add_feature(feature)

步骤 7: 显示图像

最后，我们使用 plt.tight_layout() 来自动调整布局，并显示最终结果。

# 显示图像plt.tight_layout()plt.show()

结果展示

通过上述代码，我们生成了一个展示 「南海区域」 地形的地图，并在其中嵌入了一个小地球视图，红框标出了主图显示区域。可以通过修改 lonmin、lonmax、latmin 和 latmax 的值来绘制其他海域的地图。