forked from diegormsouza/nwp-python-jul-2021
-
Notifications
You must be signed in to change notification settings - Fork 0
/
script_16.py
121 lines (92 loc) · 5.65 KB
/
script_16.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
#-----------------------------------------------------------------------------------------------------------
# INPE / CPTEC Training: NWP Data Processing With Python - Script 16: Barbs
# Author: Diego Souza
#-----------------------------------------------------------------------------------------------------------
import pygrib # Provides a high-level interface to the ECWMF ECCODES C library for reading GRIB files
import matplotlib.pyplot as plt # Plotting library
import cartopy, cartopy.crs as ccrs # Plot maps
import cartopy.io.shapereader as shpreader # Import shapefiles
import numpy as np # Scientific computing with Python
import matplotlib # Comprehensive library for creating static, animated, and interactive visualizations in Python
#-----------------------------------------------------------------------------------------------------------
# Select the extent [min. lon, min. lat, max. lon, max. lat]
extent = [-93.0, -60.00, -25.00, 18.00]
# Open the GRIB file
grib = pygrib.open("gfs.t00z.pgrb2full.0p50.f000")
#-----------------------------------------------------------------------------------------------------------
# Select the variable
ucomp = grib.select(name='U component of wind', typeOfLevel = 'isobaricInhPa', level = 850)[0]
# Get information from the file
init = str(ucomp.analDate) # Init date / time
run = str(ucomp.hour).zfill(2) # Run
ftime = str(ucomp.forecastTime) # Forecast hour
valid = str(ucomp.validDate) # Valid date / time
print('Init: ' + init + ' UTC')
print('Run: ' + run + 'Z')
print('Forecast: +' + ftime)
print('Valid: ' + valid + ' UTC')
# Read the data for a specific region
ucomp, lats, lons = ucomp.data(lat1=extent[1],lat2=extent[3],lon1=extent[0]+360,lon2=extent[2]+360)
#-----------------------------------------------------------------------------------------------------------
# Select the variable
vcomp = grib.select(name='V component of wind', typeOfLevel = 'isobaricInhPa', level = 850)[0]
# Read the data for a specific region
vcomp = vcomp.data(lat1=extent[1],lat2=extent[3],lon1=extent[0]+360,lon2=extent[2]+360)[0]
#-----------------------------------------------------------------------------------------------------------
# Select the variable
prmls = grib.select(name='Pressure reduced to MSL')[0]
# Read the data for a specific region
prmls = prmls.data(lat1=extent[1],lat2=extent[3],lon1=extent[0]+360,lon2=extent[2]+360)[0]
#-----------------------------------------------------------------------------------------------------------
# Calculate the wind speed
ws = np.sqrt(ucomp**2 + vcomp**2)
# Convert to hPa
prmls = prmls / 100
#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(8,8))
# Use the Cilindrical Equidistant projection in cartopy
ax = plt.axes(projection=ccrs.PlateCarree())
# Define the image extent
img_extent = [extent[0], extent[2], extent[1], extent[3]]
ax.set_extent([extent[0], extent[2], extent[1], extent[3]], ccrs.PlateCarree())
# Add a shapefile
# https://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2019/Brasil/BR/br_unidades_da_federacao.zip
shapefile = list(shpreader.Reader('BR_UF_2019.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='gray',facecolor='none', linewidth=0.3)
# Add coastlines, borders and gridlines
ax.coastlines(resolution='10m', color='black', linewidth=0.8)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black', linewidth=0.5)
gl = ax.gridlines(crs=ccrs.PlateCarree(), color='gray', alpha=1.0, linestyle='--', linewidth=0.25, xlocs=np.arange(-180, 180, 5), ylocs=np.arange(-90, 90, 5), draw_labels=True)
gl.top_labels = False
gl.right_labels = False
# Define de contour interval
data_min = 990
data_max = 1050
interval = 2
levels = np.arange(data_min,data_max,interval)
# Create a custom color palette
colors = ["#310d00", "#631b00", "#942800", "#c53500", "#fd6123", "#fb824e", "#faa679", "#f8c8a3", "#f6f6e1", "#f8f8e7", "#fafaee", "#fcfcf5", "#fefefc", "#e8eef5", "#d9e2ef", "#c9d7e9", "#bacce2", "#a9bbd9", "#95a1c9", "#8187b9", "#6d6da8", "#595398", "#463b87", "#382f6c", "#2a2351", "#1c1836", "#0e161b"]
cmap = matplotlib.colors.ListedColormap(colors)
cmap.set_over('#000000')
cmap.set_under('#000000')
# Plot the image
img1 = ax.contourf(lons, lats, prmls, cmap=cmap, levels=levels, extend='both')
img2 = ax.contour(lons, lats, prmls, colors='black', linewidths=0.1, levels=levels)
ax.clabel(img2, inline=1, inline_spacing=0, fontsize=10,fmt = '%1.0f', colors= 'black')
# Create a flag to determine which barbs are flipped
flip_flag = np.zeros((ucomp.shape[0],ucomp.shape[1]))
# All flags below the equator will be flipped
flip_flag[lats < 0] = 1
# Plot the barbs
img3 = ax.barbs(lons[::4,::4], lats[::4,::4], ucomp[::4,::4], vcomp[::4,::4], length = 5.0, sizes = dict(emptybarb=0.0, spacing=0.2, height=0.5), linewidth=0.8, pivot='middle', barbcolor='gray', flip_barb = flip_flag[::4,::4])
# Add a colorbar
plt.colorbar(img1, label='PSML (hPa)', orientation='vertical', pad=0.05, fraction=0.05)
# Add a title
plt.title('GFS: PMSL + Winds (850 hPa)' , fontweight='bold', fontsize=10, loc='left')
plt.title('Valid: ' + valid, fontsize=10, loc='right')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('image_16.png', bbox_inches='tight', pad_inches=0, dpi=100)
# Show the image
plt.show()