Learning Objectives:
Master python to draw wind pressure and wind direction overlay
Modified from the following two articles
- https://blog.csdn.net/qiuylulu/article/details/115304020
- https://huaweicloud.csdn.net/638079ebdacf622b8df88733.html?spm=1001.2101.3001.6650.4 &utm_medium=distribute.pc_relevant.none-task-blog-2~default~CTRLIST~activity-4-1198797 64-blog- 115304020.235^v36^pc_relevant_default_base3 &depth_1-utm_source=distribute.pc_relevant.none-task-blog-2~default~CTRLIST~activity-4-119879764-blog-115304020.235^v36^pc_re levant_default_base3 &utm_relevant_index=7
Drawing results
–
Corresponding drawing code
###Import library package
#matplotlib is used to make pictures, including pyplot, ticker and other functions for drawing pictures
#numpy simple calculations
#cartopy drawing a map
#netCDF4 read NC format data
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.pyplot as plt
#In addition, in order to control all the font size, axis thickness, font and other information in the picture to be consistent, so as to make the picture more beautified, the following settings are uniformly performed:
mpl.rcParams["font.family"] = 'Arial' #default font type
mpl.rcParams["mathtext.fontset"] = 'cm' # Math text font
mpl.rcParams["font.size"] = 12
mpl.rcParams["axes.linewidth"] = 1
#Read the file, here the pressure file of ERA5 is used to select the uv wind field on 850hpa, and the single file represents the parameters of the surface.
f1=nc.Dataset('C:\Users\huain\Desktop\msl.nc')#The selected data set downloaded in ERA5 ERA5 hourly data on single levels from 1940 to present, select Mean sea level pressure
f2=nc.Dataset('C:\Users\huain\Desktop\pha.nc')#The selected data set downloaded in ERA5 ERA5 hourly data on pressure levels from 1940 to present, select u, v wind speed of 850hpa
#Reading variables requires msl (sea level air pressure), lat, lon, time, u, v and other information.
#Note: If you don't know the variable name, you can print(f1) to see what variables are there.
In #ERA5, lat is arranged from north to south like 50->30.
#lat, lon are one-dimensional functions, msl is a three-dimensional function of (time, lat, lon), u850 and v850 are four-digit functions of (time, layer number, lat, lon), and only two-dimensional results are required for drawing. So the follow-up data must be selected and processed.
msl=f1.variables["msl"][:]/100 #pa->hpa
lat=f1.variables['latitude'][:]
lon=f1.variables['longitude'][:]
time=f1.variables['time'][:]
u850=f2.variables["u"][:]
v850=f2.variables["v"][:]
#First of all, to establish two-dimensional lon2d, and lat2d, the longitude and latitude information should be two-dimensional when drawing the coloring map later.
lon2d, lat2d = np. meshgrid(lon, lat)
msl_all=np.mean(msl[:,:,:],axis=0)
u850_aim=np. mean(u850[:,:,:], axis=0)
v850_aim=np.mean(v850[:,:,:], axis=0)
#Using the simplest lat, lon projection, the picture size is 10*8, so the canvas is established and the axis ax of the picture is obtained.
proj = ccrs.PlateCarree(central_longitude=130) #China is left
fig = plt.figure(figsize=(10,8),dpi=550) # create canvas
ax = fig.subplots(1, 1, subplot_kw={<!-- -->'projection': proj}) # create subplots
u_all=u850_aim
v_all=v850_aim
#----------- draw the map ------------------------------------ -------
ax.add_feature(cfeature.LAND.with_scale('50m'))####add land######
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),lw=2)#####add coastline#########
ax.add_feature(cfeature.OCEAN.with_scale('50m'))######add ocean########
#-----------Add latitude and longitude --------------------------------- ---
extent=[90,180,-30,50]##Latitude and longitude range
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0., color='k', alpha=0.5, linestyle='--')
gl.top_labels = False ##Turn off the coordinate display on the upper side
gl.right_labels = False ##Close the right coordinate display
gl.xformatter = LONGITUDE_FORMATTER ##Convert coordinate scale to latitude and longitude style
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1] + 1, 10))
gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3] + 1, 10))
gl.xlabel_style={<!-- -->'size':10}
gl.ylabel_style={<!-- -->'size':10}
ax.set_extent(extent) #Display the selected area
#First of all, levels define the range of the color scale of the coloring map, between 1004-1032, with an interval of 1 for a color.
#ax.contourf draws a coloring map of sea level air pressure. The parameters lon2d, lat2d, and msl_all are all two-dimensional arrays. The colormap chooses 'Spectral_r'. This depends on individual needs. There are many resources on the Internet.
#ax.quiver is to draw the wind field, and also need to input lon2d, lat2d, u, v information (at this time, uv is a three-dimensional array, you need to select the number of layers you need to become two-dimensional), scale is used to adjust The size of the drawn wind field, the larger the value, the smaller the wind line. width is the thickness of the wind line, headwidth and headlength are the thickness and length of the arrow. uv and lat, lon are displayed every 5 to make the wind field map less dense, which can be adjusted according to your own needs.
#The most important thing is that you must add transform=ccrs.PlateCarree(), the map in the picture needs to be transformed, so that the map and the two-dimensional map directly made from the data correspond to each other.
#Draw the title of each subgraph, pad represents the distance of the title from the axis.
#The coloring map and wind field map are returned, so the function is completed, and only need to pass parameters to call later.
levels = np.arange(1004, 1032 + 1, 1)
cb = ax.contourf(lon2d,lat2d,msl_all, levels=levels, cmap='bwr',transform=ccrs.PlateCarree())
cq =ax.quiver(lon2d[::10,::10],lat2d[::10,::10],u_all[::10,::10],v_all[::10,::10], color='k',scale=200,zorder=10,width=0.003,headwidth=3,headlength=4.5,transform=ccrs.PlateCarree())
#Arrow (q, position left and right, position up and down, wind speed to be displayed (10m/s), label, label is on the right side of the arrow, with axes as the coordinate system, font size)
ax.quiverkey(cq, 0.85, 1.02, 10, '10 m/s', labelpos='E', coordinates='axes', fontproperties={<!-- -->'size ':12})
#Because it is a coloring map, the condition must be colorbar. ax=(ax[0][0],ax[0][1],ax[1][0],ax[1][1]), direction To be placed vertically next to cb3, the range of tick display is 1004-1032, displayed every 4th. aspect is the aspect ratio of the colorbar, and shrink is the percentage to display. pad is the distance from the image.
The pad in #ax.tick_params means the distance between the number on the color scale and the color scale, and the length and width are the length and thickness of the color scale tick.
# draw colorbar
cb = fig.colorbar(cb,ax=(ax),orientation='vertical',ticks=np.arange(1004,1032 + 4,4),aspect=30,shrink=0.9,pad=0.02)
cb.set_label('Pressure', fontsize=18)
#--------------------add text--------------------------- -
ax.text(0.85,0.87,'L',color='red',fontsize=15,weight='bold',va='bottom',ha='center', transform=ax.transAxes)
ax.text(0.92,0.12,'H',color='red',fontsize=15,weight='bold',va='bottom',ha='center', transform=ax.transAxes)
#title
plt.suptitle('2023_05_12_23:00 850hPa Winds and Pressure', fontsize=14, y=0.93)
plt.savefig("Figure2.png")