forked from diegormsouza/nwp-python-jul-2021
-
Notifications
You must be signed in to change notification settings - Fork 0
/
script_04.py
68 lines (54 loc) · 3.08 KB
/
script_04.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
#----------------------------------------------------------------------------------------------------------------------
# INPE / CPTEC Training: NWP Data Processing With Python - Script 4: Adding a Map with Cartopy
# 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 numpy as np # Scientific computing with Python
#----------------------------------------------------------------------------------------------------------------------
# Open the GRIB file
grib = pygrib.open("gfs.t00z.pgrb2full.0p50.f000")
# Select the variable
grb = grib.select(name='2 metre temperature')[0]
# Get information from the file
init = str(grb.analDate) # Init date / time
run = str(grb.hour).zfill(2) # Run
ftime = str(grb.forecastTime) # Forecast hour
valid = str(grb.validDate) # Valid date / time
print('Init: ' + init + ' UTC')
print('Run: ' + run + 'Z')
print('Forecast: +' + ftime)
print('Valid: ' + valid + ' UTC')
# Select the extent [min. lon, min. lat, max. lon, max. lat]
extent = [-93.0, -60.00, -25.00, 18.00]
# Read the data for a specific region
tmtmp, lats, lons = grb.data(lat1=extent[1],lat2=extent[3],lon1=extent[0]+360,lon2=extent[2]+360)
#----------------------------------------------------------------------------------------------------------------------
# Convert from K to °C
tmtmp = tmtmp - 273.15
#----------------------------------------------------------------------------------------------------------------------
# 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 [min. lon, max. lon, min. lat, max. lat]
img_extent = [extent[0], extent[2], extent[1], extent[3]]
# Add coastlines, borders and gridlines
ax.coastlines(resolution='10m', color='blue', linewidth=1.5)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='red', linewidth=1.0)
gl = ax.gridlines(crs=ccrs.PlateCarree(), color='green', alpha=1.0, linestyle='--', linewidth=0.50, xlocs=np.arange(-180, 180, 5), ylocs=np.arange(-90, 90, 5), draw_labels=True)
gl.top_labels = False
gl.right_labels = False
# Plot the image
img = ax.imshow(tmtmp, origin='lower', extent=img_extent, vmin=-20, vmax=48, cmap='jet')
# Add a colorbar
plt.colorbar(img, label='2 m Temperature (°C)', extend='both', orientation='vertical', pad=0.05, fraction=0.05)
# Add a title
plt.title('GFS: 2 m Temperature' , fontweight='bold', fontsize=10, loc='left')
plt.title('Valid: ' + valid, fontsize=10, loc='right')
#----------------------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('image_4.png')
# Show the image
plt.show()