# coding:utf-8
# __author__ ='Y'
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 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),720*times)
lat_new = np.linspace(min(lat),max(lat),1440*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=r'D:\\sst分布图\\未处理\\sst.png'):
# 定位到具体经纬度
map = Basemap(llcrnrlon = 0.125, llcrnrlat = -89.875, urcrnrlon = 359.875, urcrnrlat = 89.875)
lon, lat = np.meshgrid(lon_new, lat_new)
# plt.figure(figsize=(9, 7.88)) # 设置画布大小
plt.figure (figsize=(14.4,7.2 )) # 设置画布大小
ax = plt.gca()
plt.style.use('classic')
# # 绘制经纬线
map.drawparallels(np.arange(-90., 90., 30.), labels=[1,0,0,0], fontsize=10) # 纬线
map.drawmeridians(np.arange(-180., 180., 30.), labels=[0,0,0,1], fontsize=10) # 经线
map.drawcoastlines ()
#basemap有海岸线调用函数了,不用再加海岸线
# map.readshapefile(r'M:\全球海岸线\GSHHS_h_L1', name='country', color='w')
# for shp in map.country:
# poly = Polygon(xy=shp, facecolor='w') # 填充
# ax.add_patch(poly)
# 添加Colorbar
# cbar = plt.get_cmap('rainbow')
colormesh = map.pcolormesh(lon, lat, sst_new)
#横着
# cbar = map.colorbar(colormesh, location='bottom', label="contour lines", pad="10%")
#竖着
cbar = map.colorbar (colormesh, location='right', pad="10%")
momo = nc.Dataset (openfilepath)
cbar.set_label (momo.variables['sst'].units)
# cbar.ax.tick_params (labelsize=16)
# cbar = plt.colorbar ()
# cbar.set_ticks (np.linspace (-3, 45))
# 去掉图片边框
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) # 让图片铺满画布
# 添加标题、单位
cbar.set_label("℃")
plt.title('Sea Surface Temperature')
plt.savefig(savefilepath)
# plt.show()
def setalpha(openfilepath = r'D:\\sst分布图\\未处理\\sst.png', savefilepath = r'D:\\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__':
openfilepath=r'M:\SST_data\sst.day.mean.2012.nc'
# 忽略警告
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'][:]
# 关闭nc文件
f.close ()
sst=sst[3,:,:]
lon_new=lon
lat_new=lat
sst_new=sst
sst[sst == -9.96921e+36] = 0
# lon_new, lat_new, sst_new = Interpolation(lon, lat, sst, times=1)
drawing(lon_new, lat_new, sst_new, savefilepath=r'D:\\sst分布图\\未处理\\2020'+'.png')
setalpha(openfilepath=r'D:\\sst分布图\\未处理\\2020'+'.png', savefilepath=r'D:\\sst分布图\\透明化\\2020'+'.png')
print('已完成{}号sst的绘制'.format(3))
结果