资讯详情

Python读取NC格式数据绘制风场和涡度图

Python读取NC格式数据绘制风场和涡度图

  • 一、知识点
  • 二、代码拆分
    • 导入包
    • 读取数据
    • 数据分割
    • 数据处理
    • 画图设置
    • 画图
    • 出图
  • 三、完整代码
在这里插入图片描述

一、知识点

读取NC分析数据中所需的数据格式,绘制

小知识点: 1.数据读取 2.清理数据 3.设置风羽图

二、代码拆分

导入包

#处理数据的包 import xarray as xr import numpy as np #画图的包 import matplotlib.pyplot as plt #地图的包 import cartopy.crs as ccrs import my_class   ##自己写的画地图py文件 

读取数据

#打开nc文件 ds = xr.open_dataset('download_201305.nc') # print(ds) #获取经度,纬度,时间,时间的作用下面讲 lat = ds.latitude#读取维度 lon = ds.longitude#读取经度 u = ds['u']#风场U分量 v = ds['v']#风场V分量 vo = ds['vo']#涡度 

注:这里读取的数据是所有的格点数据,但我们不需要那么多的绘图,所以我们分割数据。只需要一部分

数据分割

#设置经纬度范围,不需要画多余的数据 lat_range = lat[(lat>=22) & (lat<=40)] lon_range = lon[(lon>=102) & (lon<=126)] #根据时间,经纬度得到数据,U,V,涡度 u_region = u.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000') v_region = v.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
vo_region = vo.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
#涡度乘以100000,因为是10的负5次方,方便呈现效果
vo_region = vo_region*100000

这里涡度乘以了100000,是为了方便观察,让数据绘制出来的效果更好。

数据处理

#下面是数据的简单清洗,本来没有这个,直接填色,发现有很多负值,画图太难看了
#设置条件,找到涡度大于等于2的位置
e = (vo_region >= 2)
#利用np包,把所有不满足条件的位置改写成0
vo_region = np.where(e,vo_region,0)

把涡度小于2的全部置零,这样在填色的时候小于2的地方就白了。

画图设置

#设置画布15*7
fig = plt.figure(figsize=(15,7))
#添加子图,如果是一个就是111,一行一列第一个,设置投影方式,因为要画地图了
ax = fig.add_subplot(111,projection = ccrs.PlateCarree())
#设置x,y轴标签
ax.set_xticks(np.arange(102,127,2),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(22,41,2),crs=ccrs.PlateCarree())
xticks_str = ['102', '', '106', '', '110', '', '114', '', '118','', '122', '', '126°E']
ax.set_xticklabels(xticks_str,fontsize = 18)
yticks_str = ['22 ','24 ','26 ','28 ','30 ','32 ','34 ','36 ','38 ','40°N']
ax.set_yticklabels(yticks_str,fontsize = 20)
#画地图
my_class.readshapefile('bou2_4l.shp',linewidth=1,ax=ax)

画图

#画填色图
cf_rh = ax.contourf(lon_range,#lon
                    lat_range ,#lat
                    vo_region,#vo
                    #挑了4款效果相对较好的色条,可以换着看看效果
                    cmap='hot_r',#加入色条 #hot_r lGnBu OrRd YlGnBu
                    extend='both')#设置色条为两头尖
#色标的设置
cb=fig.colorbar(cf_rh,orientation='vertical',shrink=0.75)


#画风杆
ax.barbs(lon_range[::4],#纬度降尺度,0.25变成1
         lat_range[::4],#经度降尺度,0.25变成1
         u_region[::4,::4],#u降尺度,0.25变成1
         v_region[::4,::4]#v降尺度,0.25变成1
         ,barbcolor=['k'],#颜色
         linewidth=0.5, #宽度
         length=5, #长度
         barb_increments=dict(half=2, full=4, flag=20))#风羽图设置 半根2米,一根4米,三角20米



#写a
ax.text(0.05, 0.95, '(a)',transform=ax.transAxes, fontsize=11)

各种画图参数都在代码中注释好了。

出图

plt.show()

三、完整代码

#处理数据的包
import xarray as xr
import numpy as np
#画图的包
import matplotlib.pyplot as plt
#地图的包
import cartopy.crs as ccrs
import my_class

#打开nc文件
ds = xr.open_dataset('download_201305.nc')
# print(ds)
#获取经度,纬度,时间,时间的作用下面讲
lat = ds.latitude
lon = ds.longitude
time = ds.time
u = ds['u']
v = ds['v']
vo = ds['vo']
'2013-06-27T00:00:00.000000000'
#设置经纬度范围,我们并不需要画出多余的数据
lat_range = lat[(lat>=22) & (lat<=40)]
lon_range = lon[(lon>=102) & (lon<=126)]


#根据时间,经纬度得到数据,U,V,涡度
u_region = u.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
v_region = v.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
vo_region = vo.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
#涡度乘以100000,因为是10的负5次方,方便呈现效果
vo_region = vo_region*100000
#设置画布15*7
fig = plt.figure(figsize=(15,7))
#添加子图,如果是一个就是111,一行一列第一个,设置投影方式,因为要画地图了
ax = fig.add_subplot(111,projection = ccrs.PlateCarree())



#设置x,y轴标签
ax.set_xticks(np.arange(102,127,2),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(22,41,2),crs=ccrs.PlateCarree())
xticks_str = ['102', '', '106', '', '110', '', '114', '', '118','', '122', '', '126°E']
ax.set_xticklabels(xticks_str,fontsize = 18)
yticks_str = ['22 ','24 ','26 ','28 ','30 ','32 ','34 ','36 ','38 ','40°N']
ax.set_yticklabels(yticks_str,fontsize = 20)
my_class.readshapefile('bou2_4l.shp',linewidth=1,ax=ax)





#下面是数据的简单清洗,本来没有这个,直接填色,发现有很多负值,画图太难看了
#设置条件,找到涡度大于等于2的位置
e = (vo_region >= 2)
#利用np包,把所有不满足条件的位置改写成0
vo_region = np.where(e,vo_region,0)


#画填色图
cf_rh = ax.contourf(lon_range,#lon
                    lat_range ,#lat
                    vo_region,#vo
                    #用了系统自带的色条,文章中的色条与你给的数据不匹配,效果不好
                    #挑了4款效果相对较好的色条,可以换着看看效果
                    cmap='hot_r',#加入色条 #hot_r lGnBu OrRd YlGnBu
                    extend='both')#设置色条为两头尖
#色标的设置
cb=fig.colorbar(cf_rh,orientation='vertical',shrink=0.75)


#画风杆
ax.barbs(lon_range[::4],#纬度降尺度,0.25变成1
         lat_range[::4],#经度降尺度,0.25变成1
         u_region[::4,::4],#u降尺度,0.25变成1
         v_region[::4,::4]#v降尺度,0.25变成1
         ,barbcolor=['k'],#颜色
         linewidth=0.5, #宽度
         length=5, #长度
         barb_increments=dict(half=2, full=4, flag=20))#风羽图设置 半根2米,一根4米,三角20米



#写a
ax.text(0.05, 0.95, '(a)',transform=ax.transAxes, fontsize=11)

#出图
plt.show()

标签: 色标传感器wm03nct2a色标光电传感器红外线电梯光幕

锐单商城拥有海量元器件数据手册IC替代型号,打造 电子元器件IC百科大全!

锐单商城 - 一站式电子元器件采购平台