from netCDF4 import Datasetimport numpy as npimport matplotlib.pyplot as pltimport cartopy.crs as ccrsfrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERimport matplotlib.ticker as mticker# Open the NetCDF4 file (add a directory path if necessary) for reading:data = Dataset('F:/Rpython/lp28/data/MERRA2_300.tavg1_2d_slv_Nx.20100601.nc4', mode='r')# Run the following cell to see the MERRA2 metadata. This line will print attribute and variable information. From the 'variables(dimensions)' list, choose which variable(s) to read in below:print(data)# Read in variables:# longitude and latitudelons = data.variables['lon']lats = data.variables['lat']lon, lat = np.meshgrid(lons, lats)# 2-meter eastward wind m/sU2M = data.variables['U2M']# 2-meter northward wind m/sV2M = data.variables['V2M']# Replace _FillValues with NaNs:U2M_nans = U2M[:]V2M_nans = V2M[:]_FillValueU2M = U2M._FillValue_FillValueV2M = V2M._FillValueU2M_nans[U2M_nans == _FillValueU2M] = np.nanV2M_nans[V2M_nans == _FillValueV2M] = np.nan# Calculate wind speed:ws = np.sqrt(U2M_nans**2+V2M_nans**2)# Calculate wind direction in radians:ws_direction = np.arctan2(V2M_nans,U2M_nans)# NOTE: the MERRA-2 file contains hourly data for 24 hours (t=24). To get the daily mean wind speed, take the average of the hourly wind speeds:ws_daily_avg = np.nanmean(ws, axis=0)# NOTE: To calculate the average wind direction correctly it is important to use the 'vector average' as atan2(,) where  and  are the daily average component vectors, rather than as mean of the individual wind vector direction angle.  This avoids a situation where averaging 1 and 359 = 180 rather than the desired 0.U2M_daily_avg = np.nanmean(U2M_nans, axis=0)V2M_daily_avg = np.nanmean(V2M_nans, axis=0)ws_daily_avg_direction = np.arctan2(V2M_daily_avg, U2M_daily_avg)#Plot Global MERRA-2 Wind Speed# Set the figure size, projection, and extentfig = plt.figure(figsize=(8,4))ax = plt.axes(projection=ccrs.Robinson())ax.set_global()ax.coastlines(resolution="110m",linewidth=1)ax.gridlines(linestyle='--',color='black')# Plot windspeed: set contour levels, then draw the filled contours and a colorbarclevs = np.arange(0,19,1)plt.contourf(lon, lat, ws_daily_avg, clevs, transform=ccrs.PlateCarree(),cmap=plt.cm.jet)plt.title('MERRA-2 Daily Average 2-meter Wind Speed, 1 June 2010', size=14)cb = plt.colorbar(ax=ax, orientation="vertical", pad=0.02, aspect=16, shrink=0.8)cb.set_label('m/s',size=12,rotation=0,labelpad=15)cb.ax.tick_params(labelsize=10)plt.savefig('F:/Rpython/lp28/plot29.png',dpi=1200)plt.show()

# The filled contours show the wind speed. The "quiver" function is used to overlay arrows to show the wind direction. The length of the arrows is determined by the wind speed.# Set the figure size, projection, and extentfig = plt.figure(figsize=(9,5))ax = plt.axes(projection=ccrs.PlateCarree())ax.set_extent([-62,-38,35,54])ax.coastlines(resolution="50m",linewidth=1)# Add gridlinesgl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,linewidth=1, color='black', linestyle='--')gl.xlabels_top = Falsegl.ylabels_right = Falsegl.xlines = Truegl.xlocator = mticker.FixedLocator([-65,-60,-50,-40,-30])gl.ylocator = mticker.FixedLocator([30,40,50,60])gl.xformatter = LONGITUDE_FORMATTERgl.yformatter = LATITUDE_FORMATTERgl.xlabel_style = {'size':10, 'color':'black'}gl.ylabel_style = {'size':10, 'color':'black'}# Plot windspeedclevs = np.arange(0,14.5,1)plt.contourf(lon, lat, ws[0,:,:], clevs, transform=ccrs.PlateCarree(),cmap=plt.cm.jet)plt.title('MERRA-2 2m Wind Speed and Direction, 00Z 1 June 2010', size=16)cb = plt.colorbar(ax=ax, orientation="vertical", pad=0.02, aspect=16, shrink=0.8)cb.set_label('m/s',size=14,rotation=0,labelpad=15)cb.ax.tick_params(labelsize=10)# Overlay wind vectorsqv = plt.quiver(lon, lat, U2M_nans[0,:,:], V2M_nans[0,:,:], scale=420, color='k')plt.savefig('F:/Rpython/lp28/plot29.1.png',dpi=1200)plt.show()