如图所示，使用python通过循环读取日资料绘制区域气候态年平均的海表流场二维分布图，并用流速进行填色
使用的库如下：
import xarray as xr import numpy as np import matplotlib.pyplot as plt import glob import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter import pandas as pd import calendar import os
1、对文件进行循环读取
现有2020年一年366个日数据，以月份为命名放在12个文件夹中，每个文件夹内是各个月内的日数据。
其中，每个文件的命名是有规律可循的，放在代码里可以理解为：
一段相同的字符串+年份+月份+日期+一段相同的字符串
这有利于我们写for循环对每个文件内的每个数据进行读取处理，代码如下：
for year in range(2020,2021): for month in range(1,13): for day in range(1,calendar.monthrange(year,month)[1]+1): fn='%02.0f'%month+\ '\\'+'dt_global_allsat_phy_l4_'+'%04.0f'%year+'%02.0f'%month+'%02.0f'%day+'_20210726.nc'
2、对数据做气候态年平均处理
首先需要理清的思路是，我们在上一步写好循环后，每次都是读取一个文件，而如何对数据做年平均呢？
- 试想一下，如果我们有一个三维数据：时间x经度x纬度，其中时间的大小（size）是366，那么做年平均很简单，只需要对时间维进行平均即可。
- 基于这样的逻辑，我们知道循环中每次读取的数据是一个三维的矩阵（时间x经度x纬度：time x longitude x latitude），但是也可以看作是二维的（因为时间维只有一个数据:1 x longitude x latitude），所以我们只需要提前准备好一个三维的全部是nan值的矩阵，通过循环将每个日数据存到我们构建的三维矩阵中，之后再对时间维做平均即可。
步骤如下：
- 1、构建三维矩阵
- 2、写循环读取数据，保存数据到构建的矩阵
- 3、对保存的三维矩阵的时间维做平均
tips：
- 写个小提醒，由于不同数据的分辨率不同，不同设备的内存大小不同，而提前构建一个三维矩阵需要消耗系统的内存，越细的网格、越高的分辨率占用的内存越大，大家可以先选取一个小区域，测试熟悉一下代码，再根据需求自行构建矩阵。
我这里选的是黑潮区域，经纬度大小为：312x300
代码如下：
sla_year=np.zeros((366,312,300)) sla_year[:]=np.nan i=0 for year in range(2020,2021): for month in range(1,13): for day in range(1,calendar.monthrange(year,month)[1]+1): fn='%02.0f'%month+\ '\\'+'dt_global_allsat_phy_l4_'+'%04.0f'%year+'%02.0f'%month+'%02.0f'%day+'_20210726.nc' if(os.path.exists(path+fn)): filelist=glob.glob(path+fn) file1=filelist[0] data=xr.open_dataset(file1).sel(latitude=slice(-13,65),longitude=slice(100,175)) sla=data.sla sla_year[i,:,:]=sla i=i+1 sla_mean=np.nanmean(sla_year,axis=0)
3、绘制海表流场+填色
绘制海表流场，其实可以通过流速的两个分量，绘制箭头啦！可以通过
函数进行绘制，大家可以根据数据的分辨率、选取区域等设置箭头的间隔、粗细、长短、字体的大小。
全部代码如下：
import xarray as xr import numpy as np import matplotlib.pyplot as plt import glob import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter import pandas as pd import calendar import os path=r'E:/BaiduNetdiskDownload/2020/' sla_year=np.zeros((366,312,300)) sla_year[:]=np.nan u_year=np.zeros((366,312,300)) u_year[:]=np.nan v_year=np.zeros((366,312,300)) v_year[:]=np.nan #### ====================================================================================== #### ==================================read data ========================================== #### ====================================================================================== i=0 for year in range(2020,2021): for month in range(1,13): for day in range(1,calendar.monthrange(year,month)[1]+1): # for day in range(1,4): #dt_global_allsat_phy_l4_20200101_20210726 fn='%02.0f'%month+\ '\\'+'dt_global_allsat_phy_l4_'+'%04.0f'%year+'%02.0f'%month+'%02.0f'%day+'_20210726.nc' if(os.path.exists(path+fn)): filelist=glob.glob(path+fn) file1=filelist[0] data=xr.open_dataset(file1).sel(latitude=slice(-13,65),longitude=slice(100,175)) x=data.longitude.data y=data.latitude.data time=data.time lons,lats=np.meshgrid(x,y) sla=data.sla u=data.ugos.data v=data.vgos.data sla_year[i,:,:]=sla u_year[i,:,:]=u v_year[i,:,:]=v i=i+1 #### cal year-mean sla_mean=np.nanmean(sla_year,axis=0) u_mean=np.nanmean(u_year,axis=0) v_mean=np.nanmean(v_year,axis=0) w=np.sqrt(u_mean*u_mean+v_mean*v_mean) #### ====================================================================================== #### ===================================draw ============================================== #### ====================================================================================== fig = plt.figure(figsize=(10, 9),dpi=150) ax=fig.add_subplot(111,projection=ccrs.PlateCarree()) ax.coastlines() plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', \ edgecolor='grey', facecolor='white',zorder=2)) ax.set_xticks(np.arange(100, 175, 30),crs=ccrs.PlateCarree()) ax.set_yticks(np.arange(-10, 65, 15),crs=ccrs.PlateCarree()) ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label =False))#经度0不加标识 ax.yaxis.set_major_formatter(LatitudeFormatter()) # # # contourf \ quiver \colorbar cmap=plt.get_cmap('RdYlBu_r')#'RdYlBu_r' cb=ax.contourf(x,y,w,cmap=cmap,levels=np.arange(0,1,0.01)) ax.set_extent([100,175,-10.125,65]) ax.quiver(lons[::6,::6], lats[::6,::6], u_mean[::6,::6]*1.5,v_mean[::6,::6]*1.5,pivot='mid',\ headwidth=2,width=0.002,scale=25,transform=ccrs.PlateCarree(),color='k',angles='xy',zorder=1) # # # label \labelsize ax.tick_params(which='major', direction='out', length=10, width=0.99, pad=0.2, bottom=True, left=True, right=False, top=False) ax.set_title('Sea surface flow field_'+pd.to_datetime(time).strftime("%Y_").values[0],pad=10,fontsize=20) ax.set_xlabel('Longtitude',fontsize=20) ax.set_ylabel('Latitude',fontsize=20) plt.colorbar(cb) plt.show()
一个努力学习python的海洋人 水平有限，欢迎指正！！！ 欢迎评论、收藏、点赞、转发、关注。 关注我不后悔，记录学习进步的过程~~
