相关文章推荐
  1. 读取全球地形数据
  2. 绘制全球地形图

ETOPO2v2是美国国家海洋和大气管理局(NOAA)下属的国家地球物理数据中心(NGDC)开发的全球地形模型,包括全球陆地和海洋的地形,分辨率为2分。

  1. 读取全球地形数据

数据格式为netCDF(.nc),可以考虑使用xarray库来读取nc文件。

import xarray as xr    # 导入xarray库
ds = xr.open_dataset('ETOPO2v2c_f4.nc')    # 读取全球地形数据
ds   # 显示nc文件信息
  • 读取ETOPO2v2c_f4.nc文件信息可知,x表示经度,范围为-180~180度;y表示纬度,范围为-90~90度;z表示海拔。
  1. 绘制全球地形图

地形图的绘制使用Basemap包。Basemap是Python可视化库Matplotlib下的一个工具包,主要功能是绘制二维地图,是常用的地理数据可视化工具。尽管Basemap逐渐被Cartopy所取代,但个人认为某些地方Basemap使用起来比Cartopy更加方便好用。

  • 导入绘图所需的库
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
  • 准备绘图数据
# 准备用于绘图的数据
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))  # 纬度
lon, lat = np.meshgrid(lon, lat)  # 构建经纬网
dem = ds['z'].data  # DEM数据
    # 设置地图全局属性
    plt.rcParams['font.sans-serif'] = ['Times New Roman']  # 设置整体的字体为Times New Roman
    plt.figure(figsize=(10, 6), dpi=600)  # 设置大小和分辨率
    # 创建底图,设置地图投影为World Plate Carrée,分辨率为高分辨率,地图范围为全球
    m = Basemap(projection='cyl', resolution='h', llcrnrlon=-180, llcrnrlat=-90, urcrnrlon=180, urcrnrlat=90)
    # 设置地图经纬线,并只在左端和底端显示
    m.drawmeridians(np.arange(-180, 181, 30), labels=[0, 0, 0, 1], fontsize=10, linewidth=0.8, color='silver')  # 经线
    m.drawparallels(np.arange(-90, 91, 30), labels=[1, 0, 0, 0], fontsize=10, linewidth=0.8, color='silver')  # 纬线
    # 绘制地图
    levels = [-8000, -6000, -4000, -2000, -1000, -200, -50, 0, 50, 200, 500, 1000, 1500, 2000, 3000, 4000, 
              5000, 6000, 7000, 8000]  # 创建分级
    color = ['#084594', '#2171b5', '#4292c6', '#6baed6', '#9ecae1', '#c6dbef', '#deebf7', '#006837', 
             '#31a354', '#78c679', '#addd8e', '#d9f0a3', '#f7fcb9', '#c9bc87', '#a69165', '#856b49', 
             '#664830', '#ad9591', '#d7ccca']  # 设置色带
    m.contourf(lon, lat, dem, levels=levels, extend='both', colors=color)  # 绘图,并设置图例两端显示尖端
    # 设置图例
    cb = m.colorbar(location='bottom', pad=0.35)  # 图例在底端显示
    cb.set_ticks(levels)  # 设置色带刻度
    cb.ax.tick_params(labelsize=10)  # 刻度字号大小
    cb.set_label('Global Elevation (meter)', fontsize=12)  # 设置图例名称和字体大小
    # 保存图片并显示
    plt.savefig('global_elevation.jpg', dpi=600, bbox_inches='tight', pad_inches=0.1)  # 输出地图,并设置边框空白紧密
    plt.show()  # 显示地图
# -*- encoding: utf-8 -*-
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
# 主函数
if __name__ == '__main__':
    # 读取全球地形数据
    ds = xr.open_dataset('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))  # 纬度
    lon, lat = np.meshgrid(lon, lat)  # 构建经纬网
    dem = ds['z'].data  # DEM数据
    # 设置地图全局属性
    plt.rcParams['font.sans-serif'] = ['Times New Roman']  # 设置整体的字体为Times New Roman
    plt.figure(figsize=(10, 6), dpi=600)  # 设置大小和分辨率
    # 创建底图,设置地图投影为World Plate Carrée,分辨率为高分辨率,地图范围为全球
    m = Basemap(projection='cyl', resolution='h', llcrnrlon=-180, llcrnrlat=-90, urcrnrlon=180, urcrnrlat=90)
    # 设置地图经纬线,并只在左端和底端显示
    m.drawmeridians(np.arange(-180, 181, 30), labels=[0, 0, 0, 1], fontsize=10, linewidth=0.8, color='silver')  # 经线
    m.drawparallels(np.arange(-90, 91, 30), labels=[1, 0, 0, 0], fontsize=10, linewidth=0.8, color='silver')  # 纬线
    # 绘制地图
    levels = [-8000, -6000, -4000, -2000, -1000, -200, -50, 0, 50, 200, 500, 1000, 1500, 2000, 3000, 4000, 
              5000, 6000, 7000, 8000]  # 创建分级
    color = ['#084594', '#2171b5', '#4292c6', '#6baed6', '#9ecae1', '#c6dbef', '#deebf7', '#006837', 
             '#31a354', '#78c679', '#addd8e', '#d9f0a3', '#f7fcb9', '#c9bc87', '#a69165', '#856b49', 
             '#664830', '#ad9591', '#d7ccca']  # 设置色带
    m.contourf(lon, lat, dem, levels=levels, extend='both', colors=color)  # 绘图,并设置图例两端显示尖端
    # 设置图例
    cb = m.colorbar(location='bottom', pad=0.35)  # 图例在底端显示
    cb.set_ticks(levels)  # 设置色带刻度
    cb.ax.tick_params(labelsize=10)  # 刻度字号大小
    cb.set_label('Global Elevation (meter)', fontsize=12)  # 设置图例名称和字体大小
    # 保存图片并显示
    plt.savefig('global_elevation.jpg', dpi=600, bbox_inches='tight', pad_inches=0.1)  # 输出地图,并设置边框空白紧密
    plt.show()  # 显示地图
  • 内容仅供大家学习参考,若有不足之处,敬请大家批评指正!
  1. xarray官方文档:http://xarray.pydata.org/en/stable/index.html
  2. Basemap官方文档:https://basemaptutorial.readthedocs.io/en/latest/index.html
  3. Matplotlib官方文档https://matplotlib.org/3.2.0/index.html
  4. 色带设置:https://colorbrewer2.org/
这里写自定义目录标题欢迎使用Markdown编辑器新的改变功能快捷键合理的创建标题,有助于目录的生成如何改变文本的样式插入链接与图片如何插入一段漂亮的代码片生成一个适合你的列表创建一个表格设定内容居中、居左、居右SmartyPants创建一个自定义列表如何创建一个注脚注释也是必不可少的KaTeX数学公式新的甘特图功能,丰富你的文章UML 图表FLowchart流程图导出与导入导出导入欢迎使用Markdown编辑器你好! 这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Mar 也有网友表示:写好PPT和做好PPT在职场上就是一种能力,一份好的PPT是内容好加视觉美观。 在平时的科研过程中,我们经常会输出一些二维的平面图,二维平面图反映某个变量在二维场景下的分布情况。 在我们常规的PPT展示中,二维平面图经常会出现在内容中,但是单单一张二维图,往往会显得比较单调。要写好PPT,首先还是要有好的内容,如何有效并酷炫的展示我们的平面二维结果?先尝试让它 %matplotlib inline import os # 下面这行代码是解决在 jupyter notebook 中导入basemap出现KerErrr:"PROJ_LIB"错误的方法。可以网上搜其他方法 os.environ['PROJ_LIB'] = r'D:\Program_Software\Anaconda\pkgs\proj4-5.2.0-ha925a31_1\L...
如果不使用`mpl_toolkits.basemap`,可以使用`cartopy`库来绘制非洲地图。下面是一个简单的示例代码,可以帮助你开始: ```python import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt # 创建绘图对象和地图投影 fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) # 设置地图范围和中心点 ax.set_extent([-25, 60, -45, 45], crs=ccrs.PlateCarree()) # 添加地图特征 ax.add_feature(cfeature.LAND) ax.add_feature(cfeature.OCEAN) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=':') ax.add_feature(cfeature.LAKES, alpha=0.5) ax.add_feature(cfeature.RIVERS) # 显示地图 plt.show() 这段代码将创建一个简单的地图,其中包括非洲大陆、海洋、海岸线、边界线、湖泊和河流等特征。你可以根据需要添加或删除特征,并使用`set_extent`方法来控制地图的范围和中心点。
 
推荐文章