nc文件曲面坐标转等距离坐标并对NC文件求四季平均值（内涵极地投影坐标系（极射赤道投影））

一、Python库介绍

1. 本次测试数据为PIOMAS逐月海冰厚度数据，数据链接如下：PIOMAS Variables on Model Grid
2. 我们以海冰厚度为例，会发现他的坐标系为一个二维矩阵

接着我们使用python 将longitude输出一下

当我们看到这种逐行和逐列均不相等的就说明这个一个曲面坐标而非一个等间距坐标，如果我们需要将多组数据混合运算尤其是曲面坐标系与等间距坐标运算，会及其麻烦，那么我今天在这里就向大家介绍一种可以将曲面坐标系化为等间距坐标系的方法。

3. 核心库函数：xesmf
   xesmf 这个库只在Linux和MAC上有，额好吧，我并不会有什么MAC那种高级的电脑，今天就讲讲怎么使用Linux安装这个库。
   首先你要有一个Linux的子系统或者虚拟系统之类的，windows子系统的安装方法在我之前的博文里面有介绍过需要的小伙伴请看这里，安装好WSL之后，我们需要准备一个Anaconda大礼包，Anacodna我们需要下载X86架构的。

Anaconda安装需要使用到Linux中的Bash命令，在这里我需要郑重说明一下，比如windows中路径为：

C:\Users\Administrator\Downloads\Anaconda3-2021.11-Linux-x86_64.sh

那么在WSL子系统中的对应路径就是：

/mnt/c/Users/Administrator/Downloads/Anaconda3-2021.11-Linux-x86_64.sh

Linux中的/mnt/c/等价于windows中的C:\且Linux中不能使用\作为各级目录分割

 bash /mnt/c/Users/Administrator/Downloads/Anaconda3-2021.11-Linux-x86_64.sh

接下来一路按照他的提示走就完事了！！

接着输入

conda install -c conda-forge xesmf

安装好之后，可能还需要安装xarray,cartopy

conda install xarray
conda install cartopy

二、代码剖析部分

1. 导入相关库

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import xesmf as xe
import cmaps

2. 读取文件

ds = xr.open_dataset('/mnt/d/CSDN测试数据集/heff.nc')
# print(ds.keys())
print(ds.longitude.data)

3. 对原始数据进行插值

ds_out = xr.Dataset(coords = {"lat": np.arange(-90, 90, 1.0), "lon": ("lon", np.arange(0,360, 1))}) # 定义输出坐标系信息，我这里为随测试的，分辨率是1°的
regridder = xe.Regridder(ds_in = ds, ds_out = ds_out, method = "bilinear") # 参数ds_in要求输入原始的Dataset，参数ds_out是输出坐标系信息， 参数method 为插值方法，我这里使用的是线性插值方法
ds = regridder(ds) # ok 这个就是插值完成了

4. 按照季节对数据求平均值并获取对应数据

ds = ds.sel(lat=slice(66,90)).groupby('time.season').mean() # 首先提取北纬66至90至今的数据，接着将数据以季节进行分割
lat1  = np.array(ds['lat'].sel(lat=slice(66,90)).values)
lon1  = np.array(ds['lon'].values)
DJF = ds.sel(season='DJF')['sivol'].values
JJA = ds.sel(season='JJA')['sivol'].values
MAM = ds.sel(season='MAM')['sivol'].values
SON = ds.sel(season='SON')['sivol'].values

三、绘制极射赤道投影图

由于我们出的是四季的图，也就是说有四张，那么我就先将一些重复的代码封装成一个函数

def maxax(ax):
    import matplotlib.path as mpath #导入matplotlib.path模块
    ax.gridlines() #添加网格线
    leftlon, rightlon, lowerlat, upperlat = (-180,180,60,90) #设置经纬度范围
    img_extent = [leftlon, rightlon, lowerlat, upperlat]
    ax.set_extent(img_extent, ccrs.PlateCarree()) #设置经纬度范围和投影坐标
    ax.add_feature(cfeature.COASTLINE.with_scale('50m')) #添加海岸线
    theta = np.linspace(0, 2*np.pi, 100) #将图片设置成圆形
    center, radius = [0.5, 0.5], 0.5 ##将图片设置成圆形
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T ##将图片设置成圆形
    circle = mpath.Path(verts * radius + center) ##将图片设置成圆形
    ax.set_boundary(circle, transform=ax.transAxes) #将图片设置成圆形

接着就是出图的核心代码

fig = plt.figure(figsize=(30,12))
ax1 = plt.subplot(2, 2, 1, projection=ccrs.NorthPolarStereo()) #设置子图
ax2 = plt.subplot(2, 2, 2, projection=ccrs.NorthPolarStereo()) #设置子图
ax3 = plt.subplot(2, 2, 3, projection=ccrs.NorthPolarStereo()) #设置子图
ax4 = plt.subplot(2, 2, 4, projection=ccrs.NorthPolarStereo()) #设置子图
maxax(ax1)
maxax(ax2)
maxax(ax3)
maxax(ax4)
ax1.set_title('DJF',fontsize=20) #给每个子图添加标题
ax2.set_title('MAM',fontsize=20) #给每个子图添加标题
ax3.set_title('JJA',fontsize=20) #给每个子图添加标题
ax4.set_title('SON',fontsize=20) #给每个子图添加标题
level = np.arange(-10,10.1,1) #色带显示范围为
c7 = ax1.contourf(lon,lat,BCC_DJF_s,levels = level, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r) 
ax2.contourf(lon,lat,BCC_MAM_s,levels = level, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
ax3.contourf(lon,lat,BCC_JJA_s,levels = level, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
ax4.contourf(lon,lat,BCC_SON_s,levels = level,extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)

# 下面的是设置colorbar位置，以及将colorbar设置为横向
l = 0.22
b = 0.005
w = 0.6
h = 0.02
rect = [l,b,w,h] 
cbar_ax = fig.add_axes(rect) 
fig.subplots_adjust(bottom=0.2)     
c=plt.colorbar(c7, cax=cbar_ax,orientation='horizontal', aspect=20, pad=0.1)
c.ax.tick_params(labelsize=14)

四、 完整代码

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import xesmf as xe
import cmaps

def maxax(ax):
    import matplotlib.path as mpath
    ax.gridlines()
    leftlon, rightlon, lowerlat, upperlat = (-180,180,60,90)
    img_extent = [leftlon, rightlon, lowerlat, upperlat]
    ax.set_extent(img_extent, ccrs.PlateCarree())
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    theta = np.linspace(0, 2*np.pi, 100)
    center, radius = [0.5, 0.5], 0.5
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
    circle = mpath.Path(verts * radius + center)
    ax.set_boundary(circle, transform=ax.transAxes)

ds = xr.open_dataset('/mnt/d/heff.H1978-2014.nc')
# print(ds.keys())
print(ds.longitude.data)

ds_out = xr.Dataset(coords = {"lat": np.arange(-90, 90, 1.0), "lon": ("lon", np.arange(0,360, 1))}) 
regridder = xe.Regridder(ds_in = ds, ds_out = ds_out, method = "bilinear")
ds = regridder(ds)

ds = ds.sel(lat=slice(66,90)).groupby('time.season').mean() 
lat  = np.array(ds['lat'].sel(lat=slice(66,90)).values)
lon  = np.array(ds['lon'].values)
DJF = ds.sel(season='DJF')['h_ice'].values
JJA = ds.sel(season='JJA')['h_ice'].values
MAM = ds.sel(season='MAM')['h_ice'].values
SON = ds.sel(season='SON')['h_ice'].values

fig = plt.figure(figsize=(16,12))
ax1 = plt.subplot(2, 2, 1, projection=ccrs.NorthPolarStereo())
ax2 = plt.subplot(2, 2, 2, projection=ccrs.NorthPolarStereo())
ax3 = plt.subplot(2, 2, 3, projection=ccrs.NorthPolarStereo())
ax4 = plt.subplot(2, 2, 4, projection=ccrs.NorthPolarStereo())
maxax(ax1)
maxax(ax2)
maxax(ax3)
maxax(ax4)
ax1.set_title('DJF',fontsize=20)
ax2.set_title('MAM',fontsize=20)
ax3.set_title('JJA',fontsize=20)
ax4.set_title('SON',fontsize=20)
level = np.arange(-2,5.1,0.2)
c7 = ax1.contourf(lon,lat,DJF,levels = level, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
ax2.contourf(lon,lat,MAM,levels = level, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
ax3.contourf(lon,lat,JJA,levels = level, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
ax4.contourf(lon,lat,SON,levels = level,extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)

l = 0.22
b = 0.005
w = 0.6
h = 0.02
rect = [l,b,w,h] 
cbar_ax = fig.add_axes(rect) 
fig.subplots_adjust(bottom=0.2)     
c=plt.colorbar(c7, cax=cbar_ax,orientation='horizontal', aspect=20, pad=0.1)
c.ax.tick_params(labelsize=14)