接續昨天已經能讀取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)
視覺化完成如下,看起來就是同一筆資料,只是用不同資料格式儲存。