基于Cartopy绘制海洋表面温度动态变化图

一、实验准备

啊我们这次需要用的库有：

  • cartopoy
  • xarray
  • imageio
  • numpy
  • pandas

如果发现cartopy安装不了，可以先在网站下载依赖库。

分别是：

  • pyproj
  • Pillow
  • Shapely
  • Cartopy

我们下载好对应Python版本和OS类型的.whl文件后，在anaconda prompt中输入以下命令：

1
pip install path\your_file.whl

即可。

关于海洋遥感数据，我们可以在https://podaac.jpl.nasa.gov/上下载数据，当然，NASA的这个网站可能需要科学上网。

二、数据处理

我们下载下来的文件是.nc文件，这个类型的文件可以用xarray或者netCDF4这两个库进行读取哦。下面我们展示xarray的读取方式。

1
2
# 导入库
import xarray as xr

值得一提的是，xarray有可能跟pytorch或者tensoflow产生冲突，尤其是torchtext库。这种情况下，可能需要对现有的依赖库进行更新。

1
2
# 读取数据
data = xr.open_dataset(d)
image-20220925212334604

此时呢，数据就被我们读取进来啦。我们这次需要用到的值为SST(sea surface temperature)，因此，我们需要获取到这个维度的数据。

1
2
3
4
5
sst = data.variables["sea_surface_temperature"][:] # 通过variables属性 选择对应维度的数据

# 获取经纬度数据
lon = data.variables['lon'][:]
lat = data.variables['lat'][:]

这些数据符不符合规范，是一个值得注意的问题。我们发现，该文件的SST数据缺失了不少值，需要做一个清洗。

1
2
3
4
5
6
7
# 获取到DataFrame方便处理
sst = np.array(sst)
d = pd.DataFrame(sst[0])

print(d.describe())
# 用平均值填充
d.fillna(d.mean(), inplace=True)

三、绘制图像

好啦我们现在要用到cartopy进行绘图啦~

1
2
3
# 之前拿到的数据是xarray类型的，我们需要转换到ndarry类型
lon,lat=np.array(lon),np.array(lat)
sst=d.values

设置我们地图的投影

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
   ax = plt.axes(projection=ccrs.PlateCarree())  # 设置投影方式
   extend=[-30,-16, -60, -46] # 设置显示范围
   
   ax.set_extent(extend)  # 设置经纬度范围
   ax.stock_img()#添加地球背景
   
   ax.add_feature(cfe.OCEAN) # 添加要素
   ax.add_feature(cfe.LAND,edgecolor='black')
   ax.add_feature(cfe.LAKES,edgecolor='black')
   ax.add_feature(cfe.RIVERS)
   rivers_10m = cfe.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m')
   ax.add_feature(rivers_10m,facecolor="None",edgecolor='b')
   
    ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.2, color='k', alpha=0.5, linestyle='--') # 添加经纬网
    
   ax.coastlines() # 设置海岸线
   ax.set_title("我也不知是啥数据")

levels = np.arange(sst.min(), sst.max()+ 1, 1) # 设置颜色分辨度
cp=ax.contourf(lon,lat,sst,cmap='Spectral_r',levels=levels) # 绘制二维格网
   plt.colorbar(cp)# 设置色带
   plt.show()
image-20220925214125007

四、动图制作

这里我们就需要使用到imageio模块了。

这里简单说一下如何通过imageio制作gif动图。

首先，我们需要通过imageio.imread()API读取图片数据，将其存放在容器列表中。再通过imageio.mimsave()进行制作。

1
2
3
4
5
6
# 定义一个创建gif的方法
 def create_gif(filelist,name,dur=1.0)->None:
    IMG=[]
    for i in filelist:
        IMG.append(imageio.imread(i))
    return imageio.mimsave(name,IMG,"GIF",duration=dur)

接着获取plt的图像资源

1
2
3
4
5
plt.savefig((val:=path+"%s.png"%time),dpi=200)
lis.append(val)
    
# 释放资源
plt.close()

最终结果如下图啦，数据是随便拉的，所以对不太上

New

完整代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfe
import xarray as xr
import imageio
import numpy as np
import pandas as pd
import os

# 创建gif图像
def create_gif(filelist,name,dur=1.0)->None:
    IMG=[]
    for i in filelist:
        IMG.append(imageio.imread(i))
    return imageio.mimsave(name,IMG,"GIF",duration=dur)

def working(d,path,time,lis)->None:
    # 数据清洗
    '''

    :param data: 数据
    :param path: 保存路径
    :param time: 名字吧
    :param lis: GIF文件夹
    :return:
    '''
    data = xr.open_dataset(d)
    sst = data.variables["sea_surface_temperature"][:]

    lon = data.variables['lon'][:]
    lat = data.variables['lat'][:]


    sst = np.array(sst)
    d = pd.DataFrame(sst[0])
    # print(d.describe())
    d.fillna(d.mean(), inplace=True)

    lon,lat=np.array(lon),np.array(lat)
    sst=d.values

    ax = plt.axes(projection=ccrs.PlateCarree())  # 设置投影方式
    extend=[-30,-16, -60, -46]
    ax.set_extent(extend)  # 设置经纬度范围
    ax.stock_img()#添加地球背景
    ax.add_feature(cfe.OCEAN)
    ax.add_feature(cfe.LAND,edgecolor='black')
    ax.add_feature(cfe.LAKES,edgecolor='black')
    ax.add_feature(cfe.RIVERS)

    rivers_10m = cfe.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m')
    ax.add_feature(rivers_10m,facecolor="None",edgecolor='b')

    ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.2, color='k', alpha=0.5, linestyle='--')


    ax.coastlines()
    ax.set_title("我也不知是啥数据")

    levels = np.arange(sst.min(), sst.max()+ 1, 1) # 设置绘画精度


    cp=ax.contourf(lon,lat,sst,cmap='Spectral_r',levels=levels)
    plt.colorbar(cp)

    # ax.set_xticks(np.arange(extend[0],extend[1],5))
    # ax.set_yticks(np.arange(extend[2],extend[3],5))

    plt.savefig((val:=path+"%s.png"%time),dpi=200)
    lis.append(val)
    # 释放资源喽

    plt.close()


if __name__ == '__main__':
    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
    # plt.rcParams['axes.unicode_minus'] = False  # Solve the minus sign problems   # 现在好像不设置也没有问题了？

    path=r"C:\"

    filelist=os.listdir(path)
    Time=[]
    filedata=[]

    for i,j in enumerate(filelist):
        if j.endswith(".nc"):
            Time.append(i)
            filedata.append(j)

    IMG=[]
    savePath=path+"SavePath\\"

    if not os.path.exists(savePath):
        os.makedirs(savePath)

    for i,j in enumerate(filedata):
        working(path+j,savePath,Time[i],IMG)

    create_gif(IMG,savePath+"New.gif",1.0)