iT邦幫忙

2022 iThome 鐵人賽

DAY 23
0
AI & Data

氣象食材系列 第 23

[Day 23]進階預報資料-數值預報模式-區域預報模式(WRF-3公里)視覺化(netcdf)

  • 分享至 

  • xImage
  •  

接續昨天已經能讀取netcdf檔案,我們來畫個10米風速圖。
另外也會使用grib2資料畫風速圖。
這兩張圖應該會一樣,因為我現在的netcdf檔案是由grib2轉換來的。

import netCDF4 as nc4
import pygrib as pb

nc4 = nc4.Dataset("day22/M-A0064-21091006-048.nc4")
grb2 = pb.open("day22/M-A0064-21091006-048.grb2")

讀取檔案之後,開始算10米風速的值

#netcdf 算10米風
ws10_nc4 = (nc4["UGRD_10maboveground"][0]**2 + nc4["VGRD_10maboveground"][0]**2)**0.5
lat_nc4 = nc4["latitude"][:]
lon_nc4 = nc4["longitude"][:]

#grib2 算10米風
upname = "u-component of wind"
vpname = "v-component of wind"
u10_grb2 = grb2.select(parameterName=upname, level=10)[0]["values"]
v10_grb2 = grb2.select(parameterName=vpname, level=10)[0]["values"]
ws10_grb2 = (u10_grb2**2 + v10_grb2**2)**0.5
var1 = grb2.select(parameterName=upname, level=10)[0] 
xgrid, ygrid = var1["Nx"],var1["Ny"]
lon_grb2, lat_grb2 = var1["longitudes"].reshape(ygrid,xgrid), var1["latitudes"].reshape(ygrid,xgrid)

取值完成後,開始視覺化

import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as crs
import cartopy.io.shapereader as shpreader
from cartopy.feature import ShapelyFeature
import numpy as np
import matplotlib.colors as mcolors
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
from datetime import datetime, timedelta
from pathlib import Path
import matplotlib.font_manager as font_manager

fpath = Path("font/msjh.ttf")
msjhfont=font_manager.FontProperties(fname=fpath,size=14)
shp_tw = ShapelyFeature(shpreader.Reader("tw_shp/COUNTY_MOI_1090820.shp").geometries() , 
                        crs.PlateCarree(), facecolor='none',edgecolor='k',alpha = 0.85)
shp_ch = ShapelyFeature(shpreader.Reader("gadm41_CHN_shp/gadm41_CHN_0.shp").geometries() , 
                        crs.PlateCarree(), facecolor='none',edgecolor='k',alpha = 0.85)
shp_jp = ShapelyFeature(shpreader.Reader("gadm41_JPN_shp/gadm41_JPN_0.shp").geometries() , 
                        crs.PlateCarree(), facecolor='none',edgecolor='k',alpha = 0.85)
shp_ph = ShapelyFeature(shpreader.Reader("gadm41_PHL_shp/gadm41_PHL_0.shp").geometries() , 
                        crs.PlateCarree(), facecolor='none',edgecolor='k',alpha = 0.85)
shp_vt = ShapelyFeature(shpreader.Reader("gadm41_VNM_shp/gadm41_VNM_0.shp").geometries() , 
                        crs.PlateCarree(), facecolor='none',edgecolor='k',alpha = 0.85)

#因為要畫兩張圖,所以把個別的資訊變成list
fig, axs = plt.subplots(1,2,figsize=(16,8),subplot_kw={"projection":crs.PlateCarree()})
wslt = [ws10_nc4,ws10_grb2]
latlt = [lat_nc4,lat_grb2] 
lonlt = [lon_nc4,lon_grb2]
ttlt = ["netcdf","grib2"]
for axsflat,ws10,lat,lon,titext in zip(axs.flatten(), wslt, latlt, lonlt, ttlt):
    isw = axsflat.contourf(lon, lat, ws10, levels=range(31),cmap='plasma',transform=crs.PlateCarree())
    for shs in [shp_tw,shp_ch,shp_jp,shp_ph,shp_vt]:
        axsflat.add_feature(shs)
    axsflat.set_title(titext, fontsize=18)
    gl = axsflat.gridlines(draw_labels=True,x_inline=False, 
                           y_inline=False,linestyle="--",linewidth=0.8)
    gl.top_labels = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabel_style = {'rotation': 0, 'rotation_mode': 'anchor'}
    gl.ylabel_style = {'rotation': 0, 'rotation_mode': 'anchor'}

cbar_ax = fig.add_axes([0.25, 0.2, 0.5, 0.04]) #[x0,y0,width,height]
cbar = fig.colorbar(isw,cax = cbar_ax ,orientation="horizontal")
cbar.set_label("風速",fontproperties=msjhfont)
cbar.ax.tick_params(labelsize=13)

視覺化完成如下,看起來就是同一筆資料,只是用不同資料格式儲存。
https://ithelp.ithome.com.tw/upload/images/20221007/20150923vO04nTDJPc.png


上一篇
[Day 22]進階預報資料-數值預報模式-區域預報模式(WRF-3公里)(netcdf)
下一篇
[Day 24]進階觀測資料-雷達觀測(Rainbow5)(1/2)
系列文
氣象食材30
圖片
  直播研討會
圖片
{{ item.channelVendor }} {{ item.webinarstarted }} |
{{ formatDate(item.duration) }}
直播中

尚未有邦友留言

立即登入留言