绘制中国近海的sst分布图_sst卫星数据

---

前言

这里使用的数据为1°×1°的NOAA卫星月平均sst数据，这是下载的网址https://psl.noaa.gov/data/gridded/tables/sst.html，进去可以选择自己想要的数据。

---

提示：这里用到的库需要提前安装，具体流程可自行百度。

绘制sst分布图

思路：先提取nc文件中存储的经纬度以及sst数据，由于数据密度较低，对其进行插值。另外，由于项目对岸线的精度要求较高，下载了全球的海岸线数据进行叠加，可视要求自行选择，如对精度要求不高，Basemap库中有自带的海岸线数据。存储图片后，对陆地部分进行透明化处理，以便贴合项目。

代码如下：

import netCDF4 as nc
from matplotlib import pyplot as plt
from matplotlib.patches import Polygon
import numpy as np
from mpl_toolkits.basemap import Basemap
from scipy import interpolate
import warnings
from PIL import Image

def readnc(openfilepath = 'D:\\sst.mon.mean.nc',month=None):
    # 忽略警告
    warnings.filterwarnings("ignore")
    # 显示全部数据
    np.set_printoptions(threshold=np.inf)
    # 读取nc文件
    f = nc.Dataset(openfilepath)
    # print(dataset.variables.keys()) # 打印变量的属性值
    # 读取数据
    lat = f.variables['lat'][:]
    lon = f.variables['lon'][:]
    time = f.variables['time']
    sst = f.variables['sst'][month]
    # 关闭nc文件
    f.close()
    return lon, lat, sst

def Interpolation(lon, lat, sst, times=1):
    sst[sst>50] = 20
    # sst=np.squeeze(sst)
    func = interpolate.interp2d(lon,lat,sst,kind='linear')
    lon_new = np.linspace(min(lon),max(lon),360*times)
    lat_new = np.linspace(min(lat),max(lat),180*times)
    sst_new = func(lon_new,lat_new)#xnew, ynew是一维的，输出znew是二维的
    return lon_new, lat_new, sst_new

def drawing(lon_new, lat_new, sst_new, savefilepath='D:\\中国近海sst分布图\\未处理\\sst.png'):
    # 定位到具体经纬度
    map = Basemap(llcrnrlon = 114.5, llcrnrlat = 20.5, urcrnrlon = 130.5, urcrnrlat = 34.5)
    lon, lat = np.meshgrid(lon_new, lat_new)
    plt.figure(figsize=(9, 7.88)) # 设置画布大小
    ax = plt.gca()
    plt.style.use('classic')
    # # 绘制经纬线
    # map.drawparallels(np.arange(-90., 91., 5.), labels=[1,0,0,0], fontsize=10)  # 纬线
    # map.drawmeridians(np.arange(-180., 181., 5.), labels=[0,0,0,1], fontsize=10)  # 经线

    # m.drawmapboundary(fill_color = 'aqua')
    # map.fillcontinents(color = 'white', lake_color = 'white')
    # map.drawcoastlines()
    map.readshapefile('C:\\Users\\陌离\\Desktop\\全球海岸线shape\\GSHHS_h_L1', name='country', color='w')
    # D:\\ArcGIS10.3\\海岸线数据\\全球海岸线\\海岸线
    # C:\\Users\\陌离\\Desktop\\全球海岸线shape\\GSHHS_h_L1
    for shp in map.country:
        poly = Polygon(xy=shp, facecolor='w') # 填充
        ax.add_patch(poly)
        # 添加Colorbar
    # cmap = plt.get_cmap('rainbow')
    colormesh = map.pcolormesh(lon, lat, sst_new)
    # cb = map.colorbar(colormesh, location='bottom', label="contour lines", pad="10%")

    # 去掉图片边框
    ax.spines['right'].set_visible(False)
    ax.spines['left'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    plt.subplots_adjust(left=0, bottom=0, right=1, top=1, hspace = 0, wspace = 0) # 让图片铺满画布
    # 添加标题、单位
    # cb.set_label("℃")
    # plt.title('Sea Surface Temperature')

    plt.savefig(savefilepath)
    # plt.show()

def setalpha(openfilepath = 'D:\\中国近海sst分布图\\未处理\\sst.png', savefilepath = 'D:\\中国近海sst分布图中国近海sst分布图\\透明化\\sst_alpha.png'):
    # 将图片中的白色改为透明色
    img = Image.open(openfilepath)  # 读取照片
    img = img.convert('RGBA')    # 转换格式，确保像素包含alpha通道
    width, height = img.size     # 长度和宽度
    for i in range(0,width):     # 遍历所有长度的点
        for j in range(0,height):       # 遍历所有宽度的点
            data = img.getpixel((i,j))  # 获取一个像素
            if (data.count(255) == 4):  # RGBA都是255，改成透明色
                img.putpixel((i,j),(255,255,255,0))
    img.save(savefilepath)  # 保存图片

if __name__ == '__main__':
    for i in range(10):
        lon, lat, sst = readnc(month=i+1548)
        lon_new, lat_new, sst_new = Interpolation(lon, lat, sst, times=20)
        drawing(lon_new, lat_new, sst_new, savefilepath='D:\\中国近海sst分布图\\未处理\\2020.'+str(i+1)+'.png')
        setalpha(openfilepath='D:\\中国近海sst分布图\\未处理\\2020.'+str(i+1)+'.png', savefilepath='D:\\中国近海sst分布图\\透明化\\2020.'+str(i+1)+'.png')
        print('已完成{}月sst的绘制'.format(i+1))
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

---

以上是最终做出来的效果图，经纬度度格网、colorbar、坐标轴根据要求隐去了，可根据需要添加。

这是第一次写博客，还有很多东西不会哈哈，以后会慢慢更新的，记录下自己学习的经历。