Python基于basemap画论文级别的多子图空间场
2024-05-06 12:49:54
前言
对于地学研究中，常见的就是画带有地图的空间场，在python中常用的就是基于basemap（好像停止更新了，不太确定）和cartopy（英国气象局开发，越来越完善）进行绘制，本文基于basemap介绍一些入门：画空间场以及多子图。
安装
几年前安装basemap是真的很麻烦，现在已经很简单了，直接一句话：
conda install basemap
主要依赖的绘图库
import matplotlib as mpl from matplotlib.colors import LogNorm from mpl_toolkits.basemap import Basemap, cm, shiftgrid, addcyclic import matplotlib.pyplot as plt
可以直接使用的画图函数 （详细注释版）
这个函数可以直接使用，画每一个子图里面的空间场信息
def Plot2DForsubplots(data2D, lat, lon, DrawingNo, _reverse = True, defaultcamp = True):'''绘制 多子图空间场 ----------data2D: np.array, 2维变量数据 [lat,lon]lat: np.array, 1维, data2D相应的纬度信息lon: np.array, 1维, data2D相应的经度信息DrawingNo: str, 标号信息_reverse: 如果data2D的纬度是是从-90:90, 则需要将纬度反向, 变为90:-90-------'''t0Plot2D= time.time()print('*****************函数Plot2D开始*********************')if _reverse:data2D_plot = data2D[::-1, :] # 纬度反向 画图要90：-90else:data2D_plot = data2D# 纬度反向 画图要90：-90# ========plot================# plt.subplots(figsize=(12, 8))mpl.rcParams['font.sans-serif'] = ['Times New Roman'] # 设置matplotlib整体用Times New Romanmpl.rcParams['font.weight'] = 'bold' # 设置matplotlib整体用Times New Romanmpl.rcParams['font.size'] = 14 # 设置matplotlib整体用Times New Romanax = plt.gca()ax.yaxis.set_tick_params(labelsize=10.0)map = Basemap(projection='cyl', llcrnrlat=np.min(lat), urcrnrlat=np.max(lat), resolution='c',llcrnrlon=np.min(lon),urcrnrlon=np.max(lon)) # 所需的地图域(度)的左下角的经度 llcrnrlat The lower left corner geographical longitude# map.drawcoastlines(color='grey',linewidth=1)map.drawparallels(np.arange(-90, 91., 30.), labels=[True, False, False, True], fontsize=8,fontproperties='Times New Roman',color='grey', linewidth=0.5) # 画平行线map.drawmeridians(np.arange(0., 360, 60), labels=[True, False, False, True], fontsize=8,fontproperties='Times New Roman',color='grey', linewidth=0.5) # 在地图上画经线# map.readshapefile('/newshps/country1', 'country',# color='grey', linewidth=1)#map.fillcontinents(color='silver', zorder=1)map.etopo() # 浮雕map.drawcoastlines()# ===========自定义color===============if defaultcamp:norm = mpl.colors.Normalize(vmin=-40, vmax=40)cax = plt.imshow(data2D_plot, cmap='seismic', extent=[np.min(lon), np.max(lon), np.min(lat), np.max(lat)], norm=norm,alpha=1) # clim:色标范围 extent：画图经纬度区域 zorder：控制绘图顺序 vmin=data2D.max(), vmax=data2D.min(),else:#######==========自己设置colorbar间隔====================# mycolors = [(0.52941, 0.7, 1), (0.52941, 0.80784, 0.92157), (0.95,0.55,0.55), (0.98039, 0.50196, 0.44706), (0.99, 0.38824, 0.3),# (1, 0.28, 0.28), (1, 0.15, 0.15), (1, 0, 0)]# my_cmap = mpl.colors.ListedColormap(mycolors)# clevs = [-12,-8,-6, -4, -2, 0, 2, 4, 6, 8, 10, 12]# print(my_cmap.N)# norm = mpl.colors.BoundaryNorm(clevs, my_cmap.N)# cax = plt.imshow(data2D_plot, cmap=my_cmap,# extent=[np.min(lon), np.max(lon), np.min(lat), np.max(lat)], norm=norm,# alpha=1) # clim:色标范围 extent：画图经纬度区域 zorder：控制绘图顺序norm = mpl.colors.Normalize(vmin=-12, vmax=12)cax = plt.imshow(data2D_plot, cmap='seismic', extent=[np.min(lon), np.max(lon), np.min(lat), np.max(lat)], norm=norm, alpha=1) # clim:色标范围 extent：画图经纬度区域 zorder：控制绘图顺序# =======添加图号===============plt.text(0, 91.5, DrawingNo, fontsize=12, fontproperties='Times New Roman')#plt.title(titles, fontsize=15, fontproperties='Times New Roman')# =======设置colorbar===============# cax = plt.axes([0.35, 0.28, 0.4, 0.02]) #前两个指的是相对于坐标原点的位置，后两个指的是坐标轴的长/宽度 这个设置的colorbar位置cbar = plt.colorbar()cbar.set_label('Temperature', fontsize=12, fontproperties='Times New Roman')ax.spines['bottom'].set_linewidth(2); ###设置底部坐标轴的粗细ax.spines['left'].set_linewidth(2); ####设置左边坐标轴的粗细ax.spines['right'].set_linewidth(2); ###设置右边坐标轴的粗细ax.spines['top'].set_linewidth(2); ###设置右边坐标轴的粗细print('***********************函数Plot2D结束, 耗时: %.3f s / %.3f mins****************' % ((time.time() - t0Plot2D), (time.time() - t0Plot2D) / 60))
完整脚本：
这个脚本是调用不断调用Plot2DForsubplots来实现画多子图，如果不需要画多子图，把循环去掉就行。
import os
import time
import code
import numpy as np
import pandas as pd
import xarray as xr
#=画图需要的包
import matplotlib as mpl
from matplotlib.colors import LogNorm
from mpl_toolkits.basemap import Basemap, cm, shiftgrid, addcyclic
import matplotlib.pyplot as pltdef Plot2DForsubplots(data2D, lat, lon, DrawingNo, _reverse = True, defaultcamp = True):'''绘制 多子图空间场----------data2D: np.array, 2维变量数据 [lat,lon]lat: np.array, 1维, data2D相应的纬度信息lon: np.array, 1维, data2D相应的经度信息DrawingNo: str, 标号信息_reverse: 如果data2D的纬度是是从-90:90, 则需要将纬度反向, 变为90:-90-------'''t0Plot2D= time.time()print('*****************函数Plot2D开始*********************')if _reverse:data2D_plot = data2D[::-1, :] # 纬度反向 画图要90：-90else:data2D_plot = data2D# 纬度反向 画图要90：-90# ========plot================# plt.subplots(figsize=(12, 8))mpl.rcParams['font.sans-serif'] = ['Times New Roman'] # 设置matplotlib整体用Times New Romanmpl.rcParams['font.weight'] = 'bold' # 设置matplotlib整体用Times New Romanmpl.rcParams['font.size'] = 14 # 设置matplotlib整体用Times New Romanax = plt.gca()ax.yaxis.set_tick_params(labelsize=10.0)map = Basemap(projection='cyl', llcrnrlat=np.min(lat), urcrnrlat=np.max(lat), resolution='c',llcrnrlon=np.min(lon),urcrnrlon=np.max(lon)) # 所需的地图域(度)的左下角的经度 llcrnrlat The lower left corner geographical longitude# map.drawcoastlines(color='grey',linewidth=1)map.drawparallels(np.arange(-90, 91., 30.), labels=[True, False, False, True], fontsize=8,fontproperties='Times New Roman',color='grey', linewidth=0.5) # 画平行线map.drawmeridians(np.arange(0., 360, 60), labels=[True, False, False, True], fontsize=8,fontproperties='Times New Roman',color='grey', linewidth=0.5) # 在地图上画经线# map.readshapefile('/home/tongxuan/newshps/country1', 'country',# color='grey', linewidth=1)#map.fillcontinents(color='silver', zorder=1)map.etopo() # 浮雕map.drawcoastlines()# ===========自定义color===============if defaultcamp:norm = mpl.colors.Normalize(vmin=-40, vmax=40)cax = plt.imshow(data2D_plot, cmap='seismic', extent=[np.min(lon), np.max(lon), np.min(lat), np.max(lat)], norm=norm,alpha=1) # clim:色标范围 extent：画图经纬度区域 zorder：控制绘图顺序 vmin=data2D.max(), vmax=data2D.min(),else:#######==========自己设置colorbar间隔====================# mycolors = [(0.52941, 0.7, 1), (0.52941, 0.80784, 0.92157), (0.95,0.55,0.55), (0.98039, 0.50196, 0.44706), (0.99, 0.38824, 0.3),# (1, 0.28, 0.28), (1, 0.15, 0.15), (1, 0, 0)]# my_cmap = mpl.colors.ListedColormap(mycolors)# clevs = [-12,-8,-6, -4, -2, 0, 2, 4, 6, 8, 10, 12]# print(my_cmap.N)# norm = mpl.colors.BoundaryNorm(clevs, my_cmap.N)# cax = plt.imshow(data2D_plot, cmap=my_cmap,# extent=[np.min(lon), np.max(lon), np.min(lat), np.max(lat)], norm=norm,# alpha=1) # clim:色标范围 extent：画图经纬度区域 zorder：控制绘图顺序norm = mpl.colors.Normalize(vmin=-12, vmax=12)cax = plt.imshow(data2D_plot, cmap='seismic', extent=[np.min(lon), np.max(lon), np.min(lat), np.max(lat)], norm=norm, alpha=1) # clim:色标范围 extent：画图经纬度区域 zorder：控制绘图顺序# =======添加图号===============plt.text(0, 91.5, DrawingNo, fontsize=12, fontproperties='Times New Roman')#plt.title(titles, fontsize=15, fontproperties='Times New Roman')# =======设置colorbar===============# cax = plt.axes([0.35, 0.28, 0.4, 0.02]) #前两个指的是相对于坐标原点的位置，后两个指的是坐标轴的长/宽度 这个设置的colorbar位置cbar = plt.colorbar()cbar.set_label('Temperature', fontsize=12, fontproperties='Times New Roman')ax.spines['bottom'].set_linewidth(2); ###设置底部坐标轴的粗细ax.spines['left'].set_linewidth(2); ####设置左边坐标轴的粗细ax.spines['right'].set_linewidth(2); ###设置右边坐标轴的粗细ax.spines['top'].set_linewidth(2); ###设置右边坐标轴的粗细print('***********************函数Plot2D结束, 耗时: %.3f s / %.3f mins****************' % ((time.time() - t0Plot2D), (time.time() - t0Plot2D) / 60))def PlotTenYearAve():'''画CEMS2（第二版数据）和ERA5逐十年平均 多子图:return:'''#=加载数据Cesm2DataAnnual = np.load('/home/Cesm2DataAnnual.npy')#1950-2014 (65, 181, 359) lat:-90:90 lon:0:360Era5DataAnnual = np.load('/home/ERA5DataAnnual.npy')#1950-2014 (65, 181,360) lat:-90:90 lon:0:360LonLR = np.arange(0, 358.75, 1)LatLR = np.arange(-90, 90.5, 1)YearInterval = 10 #每十年画一张图DrawNo = ['(a)', '(b)', '(c)','(d)', '(e)', '(f)','(g)', '(h)', '(i)','(j)', '(k)', '(l)','(m)', '(n)', '(o)','(p)', '(q)', '(i)']fig = plt.figure(figsize=(25, 18))fig.subplots_adjust(wspace=0.05)for i in range((2014-1950)//10):Yeartitles = str(1950 + i * 10) + 's'Cesm2DataEach10Year = np.mean(Cesm2DataAnnual[i*YearInterval : (i+1)*YearInterval, :, :], axis=0) # (181, 359) lat:-90:90 lon:0:358Era5DataEach10Year_ = np.mean(Era5DataAnnual[i * YearInterval: (i + 1) * YearInterval, :, :], axis=0)#(181,360) lat:-90:90 lon:0:360Era5DataEach10Year = Era5DataEach10Year_[:Cesm2DataEach10Year.shape[0], :Cesm2DataEach10Year.shape[1]]# (181, 359)Cesm2_Era5Each10Year = Cesm2DataEach10Year - Era5DataEach10Yearplt.subplot((2014-1950)//10, 3, i * 3 + 1)Plot2DForsubplots(Cesm2DataEach10Year - 273.15, LatLR, LonLR, '%s CESM2 %s' % (DrawNo[i * 3 + 0],Yeartitles), _reverse=True, defaultcamp=True)plt.subplot((2014 - 1950) // 10, 3, i * 3 + 2)Plot2DForsubplots(Era5DataEach10Year - 273.15, LatLR, LonLR, '%s ERA5 %s' % (DrawNo[i * 3 + 1],Yeartitles), _reverse=True, defaultcamp=True)plt.subplot((2014 - 1950) // 10, 3, i * 3 + 3)Plot2DForsubplots(Cesm2_Era5Each10Year, LatLR, LonLR, '%s Difference %s' % (DrawNo[i * 3 + 2],Yeartitles), _reverse=True, defaultcamp=False)plt.savefig('/homeTenYearAve.jpg', dpi=600, bbox_inches='tight')if __name__ == "__main__":t0 = time.time()# ===================================================================================print('*****************程序开始*********************' )PlotTenYearAve()print('***********************程序结束, 耗时: %.3f s / %.3f mins****************' % ((time.time() - t0), (time.time() - t0) / 60))code.interact(local=locals())
效果展示
Python基于basemap画论文级别的多子图空间场相关推荐
- 【数据处理】 python 基于Basemap地理信息可视化
数据 在某些建模比赛中,经常碰见一些地理信息的数据,它们大多都是车辆在驾驶过程中随时间产生的经纬度序列,并伴随着车辆的各种参数.由于GPS信号并不是一直都很好(如过隧道或偏远地区),无论建模的问题 ...
- Python下basemap画出的各种地图
刚接触Python的basemap库时,被它所能产生的效果震撼了. 但是在深入的学习时发现网上很难找到系统的中文教程,仅能搜到一些博客文章里讲到的某些知识点,不成体系,就难以运用自如. 在网上看了看官 ...
- Python基于Django毕业设计论文选题系统设计
技术环境: PyCharm + Django2.2 + Python3.7 + mysql 系统有三个身份:学生,教师和管理员.学生在网站前端可以查询班级信息,查询教师风采,可以查询某个老师开设的论文 ...
- python简单代码画曲线图教程-用Python画论文折线图、曲线图？几个代码模板轻松搞定！...
前言 这几天在搞论文图,唉说实话抠图这种东西真能逼死人.坐在电脑前抠上一天越看越丑,最后把自己丑哭了-- 到了画折线图分析的时候,在想用哪些工具的时候.首先否决了excel,读书人的事,怎么能用exc ...
- python画折线图代码-用Python画论文折线图、曲线图？几个代码模板轻松搞定！
前言 这几天在搞论文图,唉说实话抠图这种东西真能逼死人.坐在电脑前抠上一天越看越丑,最后把自己丑哭了-- 到了画折线图分析的时候,在想用哪些工具的时候.首先否决了excel,读书人的事,怎么能用exc ...
- matlab 折线图_用Python画论文折线图、曲线图？几个代码模板轻松搞定！
前言 这几天在搞论文图,唉说实话抠图这种东西真能逼死人.坐在电脑前抠上一天越看越丑,最后把自己丑哭了-- 到了画折线图分析的时候,在想用哪些工具的时候.首先否决了excel,读书人的事,怎么能用exc ...
- 用Python画论文折线图、曲线图？几个代码模板轻松搞定！
前言 这几天在搞论文图,唉说实话抠图这种东西真能逼死人.坐在电脑前抠上一天越看越丑,最后把自己丑哭了-- 到了画折线图分析的时候,在想用哪些工具的时候.首先否决了excel,读书人的事,怎么能用exc ...
- windows下python安装basemap，画世界地图以及根据经纬度标点
环境安装 1.确定自己的python环境安装完毕且已配置好环境变量 2.安装geos:pip install geos 3.直接使用pip3 install pyproj 可能会安装错误,所以需要自己 ...
- Python 还原控制SCI论文算法系列1: 基于策略迭代的自适应最优控制器设计
Python 还原控制SCI论文算法系列1: 基于策略迭代的自适应最优控制器设计 文章目录 Python 还原控制SCI论文算法系列1: 基于策略迭代的自适应最优控制器设计 0.前言 1.研究问题的描 ...
最新文章
- c#分布式ID生成器
- c# 连续抓取页面内容
- 概率统计笔记：高斯威沙特分布
- 【转】Android 轻松实现语音识别
- linux中的加法函数,上下文管理练习（为加法函数计时）
- dalvik虚拟机简单介绍
- linux下tmux
- 广外大全国计算机,广外全国计算机等级考试考生人数再创新高
- 阿里P9工程师指定面试复习资料
- Simscape Multibody 多体动力学仿真教程（一）
- 浙大 PAT 甲级 1077 Kuchiguse
- Hash 表的时间复杂度为什么是 O(1)(面试版)
- [转贴][教学] 教你如何打飞机 ^_^
- 增量式分级判别回归树（IHDR）|翻译与笔记
- 【LintCode 题解】小米面试算法题：搜索旋转排序数组
- wdk与DDK有什么区别
- python可以开发app吗-惊呆！那些顶级App居然是用Python开发的
- 基于WSL搭建ESP8266开发环境
- tp5 接收图片_TP5框架实现上传多张图片的方法分析
- 移动端app开发-03-IOS 初级开发入门教程