气象科研绘图python代码及数据下载

# import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import numpy as np
import shapefile
import xarray as xr
from mpl_toolkits.basemap import Basemap
import pandas as pd
#设置画图字体的大小
plt.rcParams.update({'font.size':20})
#解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
#解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False
#设置画布和绘图区
fig = plt.figure(figsize=[32,18])
ax = fig.add_subplot(111)
data = pd.read_csv("./t_cityhour.csv")
alist1 = ["100001","2021/10/11 09:00:00", 0, "yes", 72, 10]
alist = ["100000","2021/10/11 09:00:00",  0, "yes", 137, 55]

data.loc[len(data)] =alist
data.loc[len(data)] =alist1

Z = data.pivot_table(index='latcen', columns='loncen', values='o3_1h').values
pd.DataFrame(Z).fillna(0,inplace=True)
X_unique = np.sort(data.latcen.unique())
Y_unique = np.sort(data.loncen.unique())
# X, Y = np.meshgrid(X_unique, Y_unique)

# 新建字典
nc_dict = {
    # nc文件的维度信息
    "dims": {"latitude": Z.shape[0], "time": 1, "longitude": Z.shape[1]},
    # nc文件的维度信息的坐标信息（lat,lon,time等）
    "coords": {
        "latitude": {
            "dims": ("latitude",),
            "data":X_unique,
        },
        "longitude": {
            "dims": ("longitude",),
            "data":Y_unique,
        },
        "time": {
            "dims": ("time",),
            "data":np.arange(np.datetime64('2013-01-01T00:00:00','s'),np.datetime64('2013-01-02T00:00:00','s'),np.timedelta64(24,'h'))
        },
    },
    # nc文件中的变量
    "data_vars": {
        "o3": {
            "dims": ("time", "latitude", "longitude"),
            "data":Z.reshape(1,Z.shape[0],Z.shape[1])
        }
    },
    # nc文件的全局属性信息
}

ds = xr.Dataset.from_dict(nc_dict)

lat = ds.latitude
lon = ds.longitude
X,Y=np.meshgrid(lat,lon)

data5 = (ds['o3'][0]) # 取第一时刻的数据 [0,::-1,:]表示第一个时次、纬度反向

ds = xr.open_dataset('gfs_0p25_2021101108.nc')
lat = ds.latitude
lon = ds.longitude
# data = (ds['UGRD_10maboveground'][0]) # 取第一时刻的数据 [0,::-1,:]表示第一个时次、纬度反向
data = (ds['UGRD_10maboveground'][0]) # 取第一时刻的数据 [0,::-1,:]表示第一个时次、纬度反向
data1 = (ds['VGRD_10maboveground'][0]) # 取第一时刻的数据 [0,::-1,:]表示第一个时次、纬度反向
data2 = (ds['PRMSL_meansealevel'][0]) # 取第一时刻的数据 [0,::-1,:]表示第一个时次、纬度反向
data3 = (ds['TCDC_entireatmosphere'][0]) # 取第一时刻的数据 [0,::-1,:]表示第一个时次、纬度反向
# data = (ds['LSMASK'])
#下面画
cbar_kwargs = {
    # 'orientation': 'horizontal',
    'label': 'μg/m3',
    'ticks': [0,160,200,300,400,800,1000]
}
levels = [
        -100000.0, 0.0, 16.0, 32.0, 48.0, 64.0, 80.0, 96.0, 112.0, 128.0, 144.0, 160.0, 164.0, 168.0, 172.0, 176.0, 180.0,
        184.0, 188.0, 192.0, 196.0, 200.0, 210.0, 220.0, 230.0, 240.0, 250.0, 260.0, 270.0, 280.0, 290.0, 300.0, 310.0,
        320.0, 330.0, 340.0, 350.0, 360.0, 370.0, 380.0, 390.0, 400.0, 440.0, 480.0, 520.0, 560.0, 600.0, 640.0, 680.0,
        720.0, 760.0, 800.0, 820.0, 840.0, 860.0, 880.0, 900.0, 920.0, 940.0, 960.0, 980.0, 1000.0, 10000.0
    ]
colors =[
    [ 0.16471, 0.58824, 0.08235, 1.0 ],
    [ 0.16471, 0.63922, 0.08235, 1.0 ],
    [ 0.16471, 0.73725, 0.08235, 1.0 ],
    [ 0.16471, 0.83137, 0.08235, 1.0 ],
    [ 0.22353, 0.88235, 0.07843, 1.0 ],
    [ 0.3451, 0.88627, 0.07059, 1.0 ],
    [ 0.39216, 0.88627, 0.06667, 1.0 ],
    [ 0.46667, 0.88627, 0.06275, 1.0 ],
    [ 0.54902, 0.8902, 0.05882, 1.0 ],
    [ 0.65098, 0.8902, 0.05098, 1.0 ],
    [ 0.7098, 0.89412, 0.04706, 1.0 ],
    [ 0.77255, 0.89412, 0.04314, 1.0 ],
    [ 0.83137, 0.89804, 0.03922, 1.0 ],
    [ 0.87451, 0.88627, 0.03529, 1.0 ],
    [ 0.8902, 0.8549, 0.03529, 1.0 ],
    [ 0.90588, 0.82353, 0.03137, 1.0 ],
    [ 0.91765, 0.79216, 0.03137, 1.0 ],
    [ 0.92941, 0.75686, 0.03137, 1.0 ],
    [ 0.93333, 0.71765, 0.03137, 1.0 ],
    [ 0.93725, 0.68235, 0.02745, 1.0 ],
    [ 0.94118, 0.64314, 0.02745, 1.0 ],
    [ 0.94118, 0.60392, 0.02745, 1.0 ],
    [ 0.9451, 0.56471, 0.02745, 1.0 ],
    [ 0.94902, 0.52549, 0.02745, 1.0 ],
    [ 0.95294, 0.48627, 0.02745, 1.0 ],
    [ 0.95294, 0.44706, 0.02745, 1.0 ],
    [ 0.95294, 0.41569, 0.02745, 1.0 ],
    [ 0.95686, 0.38039, 0.02745, 1.0 ],
    [ 0.95686, 0.3451, 0.02745, 1.0 ],
    [ 0.95686, 0.31373, 0.02745, 1.0 ],
    [ 0.95686, 0.28235, 0.02745, 1.0 ],
    [ 0.95686, 0.25098, 0.02745, 1.0 ],
    [ 0.95686, 0.21961, 0.02745, 1.0 ],
    [ 0.95686, 0.18824, 0.02745, 1.0 ],
    [ 0.95686, 0.15686, 0.02745, 1.0 ],
    [ 0.95686, 0.12157, 0.02745, 1.0 ],
    [ 0.95686, 0.0902, 0.03137, 1.0 ],
    [ 0.94118, 0.07451, 0.05882, 1.0 ],
    [ 0.92549, 0.0549, 0.0902, 1.0 ],
    [ 0.91373, 0.03529, 0.11765, 1.0 ],
    [ 0.89412, 0.01961, 0.16078, 1.0 ],
    [ 0.8549, 0.01961, 0.23137, 1.0 ],
    [ 0.81961, 0.01961, 0.30588, 1.0 ],
    [ 0.78431, 0.01961, 0.37647, 1.0 ],
    [ 0.74902, 0.01961, 0.45098, 1.0 ],
    [ 0.71373, 0.01961, 0.52157, 1.0 ],
    [ 0.67843, 0.01961, 0.59216, 1.0 ],
    [ 0.64314, 0.01961, 0.66667, 1.0 ],
    [ 0.60784, 0.01961, 0.72549, 1.0 ],
    [ 0.57255, 0.01961, 0.76863, 1.0 ],
    [ 0.53725, 0.01961, 0.81176, 1.0 ],
    [ 0.50196, 0.01961, 0.85882, 1.0 ],
    [ 0.47451, 0.01961, 0.83137, 1.0 ],
    [ 0.45098, 0.01961, 0.75686, 1.0 ],
    [ 0.42745, 0.01961, 0.68627, 1.0 ],
    [ 0.40784, 0.01961, 0.61569, 1.0 ],
    [ 0.38039, 0.01961, 0.54118, 1.0 ],
    [ 0.35686, 0.01961, 0.47059, 1.0 ],
    [ 0.33333, 0.01961, 0.4, 1.0 ],
    [ 0.3098, 0.01961, 0.32549, 1.0 ],
    [ 0.28627, 0.01961, 0.2549, 1.0 ],
    [ 0.23529, 0.01961, 0.1098, 1.0 ]
  ]

# cs = data5.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, colors = colors)
# cs = data3.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, colors = colors)
# cs = data2.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, colors = colors)
cs = data3.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, colors = colors)
CS3 = plt.contour( lon,lat, data2,colors = ('b'))

#
# #添加掩膜路径，白化外部的分部
# ##################################
#读取shp格式的地图，有三个文件分别为.dbf,.shp,.shx缺一不可
#这里的shp文件其实就是一堆点，把各个国家圈出来
#sf得到的是一个类，存放了全球的地理信息
sf = shapefile.Reader("country1")
#sf.shapeRecords()里面放了各个国家地区的信息
#循环的目的是单个拿出来，找到中国  shape_rec是单个
for shape_rec in sf.shapeRecords():
    #shape_rec.record内部存放了单一国家的很多信息，比如名字，货币等等
    #其中shape_rec.record[2]放的是国家名
    #可以print(shape_rec.record)看看
    if shape_rec.record[2] == 'China':
        codes = []#用来存放移动路径(画图动作)
        #shape_rec.shape是一个类，图形类
        #里面三个属性shapeType：图形类型  points图形边界坐标  parts图形起始索引
        #解释一下parts属性，一个国家的边界可能不是全连在一起的，会分为一块一块，那么就相当于多个图形，存在shp文件内就不连续，parts里面就放了每个区域的起始索引（下标）
        pts = shape_rec.shape.points
        #上文已经说过了parts的意义，设想一下，两块区域的的起始坐标中间夹的不就是一块区域的所有坐标吗，但是最后一块区域没有结束值，所有加了[len(pts)]，这就是最后一个点的索引。
        prt = list(shape_rec.shape.parts) + [len(pts)]
        #下面的循环主要作用是：建立一个绘图路径，利用区块起始点的索引生成一个画图动作
        for i in range(len(prt) - 1):
            codes += [Path.MOVETO]#点移动
            codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)#画线
            codes += [Path.CLOSEPOLY]#这块画完，循环结束，下一块
        clip = Path(pts, codes)#利用数据和路径生成一个画图动作
        clip = PathPatch(clip, transform=ax.transData)#再加入ax的变换

#生成一个Basemap类，用来变换坐标，画出合适尺寸的图
m = Basemap(llcrnrlon=72.0,#地图左边经度
    llcrnrlat=10.0,#地图下方纬度
    urcrnrlon=137.0,#地图右边经度
    urcrnrlat=55.0)#地图上方纬度

m.drawcoastlines()
m.drawmapboundary(fill_color=[193/255,211/255,255/255,1])
m.fillcontinents(color=[255/255,255/255,255/255,0.4], lake_color=[193/255,211/255,255/255,1])  # 画大洲，颜色填充

for contour in cs.collections:
        contour.set_clip_path(clip)

#读取图形文件，画中国行政边界

m.readshapefile('bou2_4l','China Map',linewidth=1.2)

######################至此图基本画完了，下面是一些修饰过程######################
#画纬度
parallels = np.arange(10,55,10)
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3])
#画经度
meridians = np.arange(70,135,10)
m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[1, 3])
#移除原来的坐标轴标签
plt.ylabel('')
plt.xlabel('')

#设置标题
ax.set_title(u' 测试分布图',color=[0.1,0.1,0.1,1],fontsize= 30)

######################修饰完毕，出图######################`
plt.savefig("测试分布图.png", dpi=300, bbox_inches='tight')
plt.show()