Contact: If you have any questions, please contact:

• E-mail: [email protected]

Learning Objectives: By the end of this activity, participants will:

• Learn how to install and use basic tools to start manipulating numerical weather preduction data with Python in your computers.
• Carrying out operations such as:

– Reading GRIB files.
– Checking the available fields in a NWP GRIB file.
– Create a basic plot visualize pixel values.
– Modify color palettes, read metadata, add a title and legend to the plot.
– Add maps and shapefiles.
– Plot contours and labels.
– Create custom color palettes.
– Soften the contours.
– Process multiple files.
– Creating animations.

During the virtual course, we will also learn:

• Data download using scripts.
• Multiple plots.
• Read several NWP fields, performing more elaborate calculations.
• Satellite Plots + NWP.
• METAR plots + Satellite + NWP.
• Among other operations!

Estimated duration of this pre-course activity: 2:00 h

Note: The number of Ads in the free version of wordpress is increasing considerably. To avoid seeing ads, I suggest installing Google Chrome’s Adblock extension:

https://chrome.google.com/webstore/detail/adblock-%E2%80%94-best-ad-blocker/gighmmpiobklfepjocnamgkkbiglidom

PREREQUISITES

For this exercise, we’ll need the following:

• Python 3.10 (our programming language)
• A “Package Manager” (to install libraries)
• An “Environment Manager” (to separate our projects)
• A Text Editor (to write our code)
• NWP data samples (the data to be manipulated)

For the first three items (“Python 3.10”, “Package Manager” and “Environment Manager”), the “Miniconda” tool will be sufficient. As for the text editor, there are many options available (Visual Studio Code, Spyder, PyCharm, Atom, Jupyter, etc), but for simplicity (and avoiding problems with editors), we will use “Notepad++”.

For NWP data, we will download GFS samples from the NOMADS server (NOAA). We can also obtain NWP data from other means (INPE’s FTP server for example).

EXECUTING SCRIPTS LOCALLY X EXECUTING SCRIPTS IN THE CLOUD

It is possible to run the Python scripts that we will see both locally (on your own computers) and in the cloud (on the Google Colab platform, for example). During this training, we will exemplify both! In this pre-course procedure, we will see how to install the necessary tools and run the scripts locally. During the course, we will see how to run the scripts in the cloud, without the need for local installation. This way, participants will have an overview of both possibilities.

With this preliminary information, let’s get to work!

IMPORTANT: If you already installed Miniconda, created a virtual environment and downloaded an Editor (following the procedure from the GOES-R activity), please jump to step number 4 (Downloading Weather Prediction Data in Grib Format).

INSTALLATION STEPS

1. Download and install Miniconda for Python 3.10 at the following link (approximatelly 60 MB):

https://docs.conda.io/en/latest/miniconda.html

Notes – Windows installation:

• During the installation, it is not necessary to check “Add Anaconda to my Path environment variable”.
• You may check “Register Anaconda as my default Python 3.10”

Notes – Linux installation:

• In the Linux installation, it is possible to choose the installation directory with the “-p” parameter. Choose the appropriate directory on your machine:

./Miniconda3-py310_23.1.0-1-Linux-x86_64.sh -p /my_directory/

Notes – Windows and Linux installation:

• The installation will take approximately 5 minutes.

CREATING A PYTHON VIRTUAL ENVIRONMENT AND INSTALLING LIBRARIES

2. Let’s create a Python environment called “workshop”. We will call it “workshop”, but you may give any name you want to your environment, as long as you use this name during the activation as we will see. We willinstall the following libraries and their dependencies in this environment:

• matplotlib: library for creating plots
• netcdf4: read and write NetCDF files
• cartopy: produce maps and other analysis of geospatial data
• boto3: integration with Amazon cloud services
• gdal: reprojection/manipulation of geospatial data
• scipy: array manipulation
• pandas: data manipulation and analysis
• pygrib: read and write GRIB files
• metpy: reading, visualization and calculations with meteorological data
• imageio: read and write a wide variety of image data

In order to create the environment, open the recently installed “Anaconda Prompt” as an Administrator and insert the following command:

conda create --name workshop -c conda-forge matplotlib netcdf4 cartopy boto3 gdal scipy pandas pygrib metpy imageio

Linux: It is not necessary to open the Anaconda Prompt, and the commands can be used in the terminal.

Where:

workshop: The name of your environment. It could be anything, like “satellite”, “goes”, etc…

matplotlib netcdf4 cartopy boto3 gdal scipy pandas pygrib metpy imageio: The libraries we wan’t to install in this environment. There are many other libraries available, these are the ones we will be using in the training.

Notes: 

• During the installation, type “y” and Enter  to proceed, when requested.
• NOTE: Creating the virtual environment and installing the libraries will take some time. Depending on the internet connection and your computer specs, it may take a considerable amount of time.

Finally, activate the newly created virtual environment with the following command:

activate workshop

Note: If you are using Linux, you can also use the following command:

source activate workshop

Extra: Although the commands below are not required for this procedure, they are nice to know for future development:

• Deactivate a conda environment: conda deactivate
• Viewing the list of environments: conda env list
• Viewing the list of packages (libraries, etc.) installed in an environment: conda list
• Deleting and existing environment (example with the “workshop” env): conda remove –name workshop –all

More information on managing environments may be found at this link.

We are ready to start using Python, however, we need an editor to write our code and we also need some NWP data samples.

DOWNLOADING AN EDITOR

3. Download and install Notepad++ at the following link (4 MB only):

https://notepad-plus-plus.org/download

Linux: You can install Notepadqq (the equivalent for Linux) with the following command:

# sudo apt-get install notepadqq

You may also use Gedit, or another Linux editor.

Note:

• You may use any editor you want (Gedit, Visual Studio Code, Spyder, PyCharm, Atom, Jupyter Notebooks, etc.), as long as you save the file with the “.py” extension.

DOWNLOADING NUMERICAL WEATHER PREDICTION DATA IN GRIB FORMAT

4. Create a directory on your machine for the files for this activity. In this example procedure we will use the VLAB\Python\NWP directory (at: C:\VLAB\Python\NWP\). This will be our working directory. Use whatever directory you like.

Let’s download a GFS sample in GRIB format from the NOMADS server. Access the following page:

https://nomads.ncep.noaa.gov/

Let’s download a sample with 0.50° resolution, time f+000. On the NOMADS server page, choose the following option: GFS 0.50 Degree – grib filter.

You will see a list of the last 10 days. No historical data is available on this server. Choose the day you want. In this example, we will choose the last day available during the development of this tutorial, March 20, 2023:

In the next two pages, we will choose the desired run (00 z in this example), and the subdirectory “atmos” (only option available).

When choosing the “atmos” subdirectory, we are forwarded to the page where we can select the following characteristics of the GFS file to be downloaded:

Desired file: f+000, f+003, f+006, f+009 … Up to f+384 (+16 days)
Desired levels: surface, 100 MB, 200 MB, etc. We can also select all (activating the “all” option).
Desired variables: temperature, etc. We can also select all (activating the “all” option).
Desired subregion: Longitudes and latitudes, if we wish a crop (activating the “make subregion” option).

In the example of this procedure, we will choose the following characteristics:

• File: f+000
• Levels: All (activating the “all” option)
• Variables: All (activating the “all” option)
• Subregion: Check the option “make subregion” and choose:
• left longitude: 260
• right longitude: 320
• top latitude: 40
• bottom latitude: 0
  

Download the file by clicking on “Start Download”, and put the downloaded file in your newly created work directory (for example: D:\VLAB\Python\NWP\)

Note: You may also find GFS data in your GEONETCast-Americas station, at the MARN-El Salvador directory. We are not going to use these files in this procedure (they have some particularities that we are going to mention during the couse).

At this point we already have our programming environment installed, our text editor installed and a GFS sample available. We are ready to start!

SCRIPT 1: KNOWING THE AVAILABLE VARIABLES

IMPORTANT: In this pre-course activity we will not explain in detail what each instruction does. The concepts will be reviewed and explained in more detail in the virtual course.

1. Open your text editor (in this procedure, “Notepad ++”), enter the code below and save it as “script_01.py” into your work directory (C:\VLAB\Python\NWP, in this example). This is the same directory where you saved the grib file in the previous step.

IMPORTANT: In line 12, insert the name of the grib file you downloaded in the previous step.

The image below shows the script open in the Notepad++ editor:

Basically, on line 10 we specify the file we wish to manipulate, on line 13, we opened the file, on line 16, we created a text file called “Variables.txt”, and between lines 17 and 21 we started a “for” loop, where for each detected variable, a new line in the text file is created with the description of the variable. In addition, each line is shown at the terminal. In line 22, we closed the text file.

In the “Anaconda Prompt”, with the “workshop” virtual environment activated, go to the C:\VLAB\Python\NWP\ (or your work directory):

cd C:\VLAB\Python\NWP

Run the “script_01.py” with the following command:

python script_01.py

When running script_01.py, you will see the variables available within the grib file directly on the terminal, and a file called “Variables.txt” will be created in your work directory:

Given the large amount of variables, it is not very practical to see in the terminal. Let’s check the text file.

When opening the text file “variables.txt”, we can check the available variables. Search for the variable “2 Metre Temperature”. Note that we can see, in addition to the name of the variable (important in the next script), the unit (k), the level (2 m), among other information.

Let’s plot the “2 Metre Temperature” variable in our next script!

SCRIPT 2: BASIC PLOT

Let’s create our second script and plot our first image. Create and save a file called “script_02.py”.

As we now wish to plot an image, besides the pygrib, let’s also import the Matplotlib library:

#-----------------------------------------------------------------------------------------------------------
import pygrib                        # Provides a high-level interface to the ECWMF ECCODES C library for reading GRIB files
import matplotlib.pyplot as plt      # Plotting library
#-----------------------------------------------------------------------------------------------------------

In the following instructions, in addition to opening the grib file as we did in the previous script, we selected the desired variable, inserting the same name as the list of variables:

# Select the variable
grb = grib.select(name='2 metre temperature')[0]

And define the region we want to read the data:

# Select the extent [min. lon, min. lat, max. lon, max. lat]
extent = [-100.0, 0.00, -40.00, 40.00]

And we read the variable, latitude and longitude arrays. We can assign the name we want to the variables, in this example we attribute as “tmtmp” (from “Two Metre Temperature”). Also note that we add 360 in the longitudes, as we specify the desired area with values between -180 and +180°, in the “extent” variable.

# Read the data for the selected extent
tmtmp, lats, lons = grb.data(lat1=extent[1],lat2=extent[3],lon1=extent[0]+360,lon2=extent[2]+360)

From there, we proceed to the plot configuration, where we chose its size (in the example below, 7×5 inches [width x height]) and its basic characteristics, such as the color palette to be used (“jet” in this case). We’ll see more details about this in the next scripts.

#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(7,5))
 
# Plot the image
plt.imshow(tmtmp, origin='lower', cmap='jet')
#----------------------------------------------------------------------------------------------------------- 

Finally, we save the plot as “image_2.png” (simply because it is the second script – you may choose the name you want) and showed the plot in an interactive window.

#----------------------------------------------------------------------------------------------------------- 
# Save the image
plt.savefig('image_2.png')

# Show the image
plt.show()

Below, the complete script so far:

Run the “script_02.py” script with the following command:

python script_02.py

When running the script, a new window will appear with our first python-made NWP plot:

• Checking the Pixel Values

Move the mouse pointer over the plot and you will see the pixels values of this particular field in Kelvin (remember the unit in the field description in the text file?) being shown in the window. In the example image below, the temperature in this pixel example is 286.5 K:

• Zoom on a given region

In order to zoom on a given region, just click on the magnifier icon in the upper part of the screen and select the region you want to zoom in.

To return to full view, click on the “Home” icon.

In addition to the interactive window, a PNG image called ‘image_02.png’ was saved in your work directory (C:\VLAB\Python\NWP in this procedure).

What can you quickly change in this script?

• The field to be read: See the text file with the variables some other variable of interest and try to plot it.
• The image to be read: Try to download other files from the NOMADS server (other dates and runs) and make a new plot.
• The size of the image generated: Try to generate images with other dimensions. What would you change to choose the size in centimeters instead of inches?

SCRIPT 3: READING METADATA, CONVERSION TO CELSIUS, COLOR PALLETE, ADDING LEGEND AND TITLE

To check the metadata available in a grib file, we can use the following command:

# Print the available keys
print("GRIB Keys :", grb.keys())

To make a quick test, add the above instruction to one of the previous scripts and see the metadata available at the terminal. Here’s an example:

#-----------------------------------------------------------------------------------------------------------
import pygrib                        # Provides a high-level interface to the ECWMF ECCODES C library for reading GRIB files
import matplotlib.pyplot as plt      # Plotting library
#----------------------------------------------------------------------------------------------------------- 

# Open the GRIB file
grib = pygrib.open("gfs.t00z.pgrb2full.0p50.f000")
 
# Select the variable
grb = grib.select(name='2 metre temperature')[0]

# Print the available keys
print("GRIB Keys :", grb.keys())

When you run the script above, you will see all the “key” information available to this field. Note that among them is the following information: “hour”, “forecastTime”, “analDate” and “validDate”:

Let’s check this information? Let’s create our third script. Create and save a file called “script_03.py”.

To read the init time, run, forecast and valid hour information, we can use the following commands:

# 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')

To convert from K to °C, we can use the following instruction:

# Convert from K to °C
tmtmp = tmtmp - 273.15

And to add a title to the plot, we can use the following command:

# Add a title
plt.title('GFS: 2 m Temperature' , fontweight='bold', fontsize=10, loc='left')
plt.title('Valid: ' + valid, fontsize=10, loc='right')

Make a test and also change the plot’s “colormap” to “Rainbow”:c

# Plot the image
plt.imshow(tmtmp, origin='lower', cmap='rainbow')

This should be your complete script by now:

Run the “script_03.py” script with the following command:

python script_03.py

The following image will appear, with a title, legend and other decorations:

• Changing the color pallete (cmap):

Now, on line 39, change from ‘rainbow’ colormap to ‘nipy_spectral’.

# Plot the image
plt.imshow(tmtmp, origin='lower', cmap='nipy_spectral')

Check at the following link the standard color palettes of Matplotlib (new color palettes can be created, as we will see later):

https://matplotlib.org/stable/tutorials/colors/colormaps.html

In addition to the interactive window, a PNG image called ‘image_03.png’ was saved in your work directory (C:\VLAB\Python\NWP in this procedure).

What can you quickly change in this script?

• The color pallete: Try to create plots with different colormaps from the link above.
• The title of the image: Try to change the size and text of the title.
• The appearance of the color bar of the caption: Try to change the orientation to “vertical” and the color bar size by modifying the value of the “fraction” variable.

Our plot is still very basic and not that interesting. Let’s improve it by adding maps and other features.

SCRIPT 4: ADDING MAPS WITH CARTOPY

Let’s add maps and latitude and longitude gridlines to our plot. To accomplish this, we will import two new libraries in this script, Cartopy and Numpy:

#-----------------------------------------------------------------------------------------------------------
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
#----------------------------------------------------------------------------------------------------------- 

To use Cartopy in our code, we need to specify the projection (“Platecaree”, or equidistant cylindrical in our case), and the plotregion [min. lon, max. lon, min. lat, max. lat]. Thus, the library will know how to position the map over the plot.

#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(10,10))

# 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]] 

And to insert maps and latitudes and longitudes in the plot, we can use the following commands:

# 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

Note that in the gridlines command, we can define the appearance of latitude and longitude lines, in addition to the desired interval (in this example, every 5 degrees), for both latitude and longitudes. Also, we disable the upper and right coordinate labels of the plot, setting them as “false”.

The last modifications required are in the imshow command, where we add the desired region, defined above. Note that we can also add the minimum and maximum thresholds to be used by color palettes:

# Plot the image
img = ax.imshow(tmtmp, origin='lower', extent=img_extent, vmin=-20, vmax=48, cmap='jet')

This should be your complete script by now:

Run the “script_04.py” script with the following command:

python script_04.py

The following image will be visualized, with the maps and gridlines of latitude and longitude:

• Viewing the latitudes, longitudes and pixel values:

Move the mouse pointer over the plot and, in addition to the pixel values of the field in ° C, you will see the pixel coordinates:

In addition to the interactive window, a PNG image called ‘image_04.png’ was saved in your work directory (C:\VLAB\Python\NWP in this procedure).

What can you quickly change in this script?

• The colors of the maps (in the variables “color” y “edgecolor”).
• The thickness of the lines (in the “linewidth” variable).
• The interval of the latitude and longitude gridlines.
• The minimum and maximum thresholds to be used by the color palette.

In this link, please check this color list with a name for each color, which can be specified in the code (“color” and “edgecolor” parameters).

Cartopy provides some custom maps. But sometimes, the maps we want to add are not available by default. That’s when adding shapefiles could be very useful. Let’s see how to do it.

SCRIPT 5: ADDING SHAPEFILES

Let’s create our fifth script. Create and save a file called “script_05.py”.

Cartopy’s standard maps are somewhat limited. We can add shapefiles, increasing the possibilities of our plots.

Please download a sample shapefile, with the world states and provinces at this link. Extract it (unzip) in your C:\VLAB\Python\NWP\ directory, where you have your scripts and sample imagery.

Following all the steps of the procedure so far, this should be the content of your work directory for now:

To add a shapefile, we need to import the “shapereader” module from the Cartopy library:

#-----------------------------------------------------------------------------------------------------------
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
#----------------------------------------------------------------------------------------------------------- 

And add shapefile to to plot with the following commands:

# Add a shapefile
shapefile = list(shpreader.Reader('ne_10m_admin_1_states_provinces.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='red',facecolor='none', linewidth=0.3)

Note that in the command “shreader.Reader” we must insert the same name of the file that has the “.shp” extension, which we have just unzipped. And in the next line, we can choose the color, thickness of the line among other characteristics.

This should be your complete script by now:

Run the “script_05.py” script with the following command:

python script_05.py

Important: This particular shapefile has a considerable size. The plot may take a while to complete. 

You should see the following plot, with the shapefiles with the states and provinces:

What can you quickly change in this script?

• The shapefile to be read. Try adding other shapefiles of your choice (for example: rivers or cities in your region).
• Shapefile characteristics. Try modifying the line colors and line thicknesses.

SCRIPT 6: PLOTTING CONTOURS AND LABELS

Let’s create our sixth script. Create and save a file called “script_06.py”.

In this script we will learn how to plot the NWP field with contours, with filling, without filling, as well as the “labels”, with the indications of the values of the contours.

The first step is to define the desired contour interval. To make it easier, we’ll use the numpy.arange function, as we use it to create the intervals of latitude and longitude lines, defining the minimum value, the maximum value and the interval:

# Define de contour interval
data_min = -20
data_max = 48 
interval = 2
levels = np.arange(data_min,data_max,interval)

to plot the contours in the defined interval, we will use the commands contourf (filled contour), contour (unfilled contour). And to show labels we will use the clabel command.

# Plot the contours
img1 = ax.contourf(lons, lats, tmtmp, cmap='jet', levels=levels, extend='both')    
img2 = ax.contour(lons, lats, tmtmp, colors='white', linewidths=0.3, levels=levels)
ax.clabel(img2, inline=1, inline_spacing=0, fontsize='10',fmt = '%1.0f', colors= 'black')

In the contourf function, we pass the latitudes, longitudes, the variable in question, the color palette and the levels. In the contour function, in this case, we pass only a color and the font size.

This should be your complete script by now:

Run the “script_06.py” script with the following command:

python script_06.py

The following image will be displayed, with the contours and values of the labels:

What can you quickly change in this script?

• Minimum and maximun values, and contour range: Test different thresholds and ranges.
• Contour characteristics: Test out different colormaps and colors on filled and unfilled contours.
• Label characteristics: Try modifying sizes and colors of labels.

• Plotting unfilled contours, with colormaps:

Let’s check how to plot only the contours without filling. Comment out the line with the contourf command and make the following changes to the contour and clabel command lines. Also comment out the colorbar instruction:

# Plot the contours
#img1 = ax.contourf(lons, lats, tmtmp, cmap='jet', levels=levels, extend='both')    
img2 = ax.contour(lons, lats, tmtmp, cmap='jet', linewidths=1.0, levels=levels)
ax.clabel(img2, inline=1, inline_spacing=0, fontsize='10',fmt = '%1.0f', colors= 'blue')

# Add a colorbar
#plt.colorbar(img1, label='2 m Temperature (°C)', orientation='vertical', pad=0.05, fraction=0.05)

Below, the new plot generated, with the unfilled contours:

SCRIPT 7: CREATING CUSTOM COLOR PALLETES

Let’s create our seventh script. Create and save a file called “script_07.py”.

We saw in previous scripts that Matplotlib has several default color palettes, which can be seen at the following link:

https://matplotlib.org/stable/tutorials/colors/colormaps.html

Despite being varied, they do not always meet the needs of the plot, and sometimes it is necessary to create custom color palettes. You can create color palettes with any number of colors, with any color. Let’s see an example.

Let’s reproduce the temperature color palette found on the following page:

http://wxmaps.org/pix/temp8

This palette is not available in matplotlib’s list of default colormaps. To reproduce it, we need to specify via code the representation of each of these colors, but in hexadecimal. How can we do that? Let’s access the following page:

https://imagecolorpicker.com/

In the webpage, click at the “Use Your Image” button:

In the next window, select the “Website Url” tab and paste the address of the webpage of interest. Click “OK”.

The webpage will automatically take a printscreen and you can click anywhere on the page image to get the hexadecimal value of the colors, thus being able to use these values in the scripts.

Now let’s insert the hex values into our script and create a custom color palette. In this example, we have 30 colors in the scale, plus one color above the upper limit and one color below the lower limit. The code below shows the colors already in hexadecimal, and the other necessary instructions:

# Create a custom color palette 
colors = ["#d3d2d2", "#bcbcbc", "#969696", "#1464d2", "#1e6eeb", "#2882f0", 
"#3c96f5", "#50a5f5", "#78b9fa", "#96d2fa", "#b4f0fa", "#1eb41e", "#37d23c", 
"#50f050", "#78f573", "#96f58c", "#b4faaa", "#c8ffbe", "#ffe878", "#ffc03c", 
"#ffa000", "#ff6000", "#ff3200", "#e11400", "#c00000", "#a50000", "#785046", 
"#8c6359", "#b48b82", "#e1beb4"]
cmap = matplotlib.colors.ListedColormap(colors)
cmap.set_over('#fadad5')
cmap.set_under('#e5e5e5')

In the “colors” list, let’s add as many colors as necessary, in hexadecimal. In the “set_over” parameter, the color of the values above the upper limit of the scale, and in “set_under” parameter, the color of the values below the lower limit of the scale.

This should be your complete script by now:

Run the “script_07.py” script with the following command:

python script_07.py

The following image will be displayed, with the custom color pallete:

• Exercise: Access the numerical weather prediction webpage form CPTEC:

https://previsaonumerica.cptec.inpe.br/

and reproduce the color scale used in the Temperatura do Ar a 2 metros (2 metre Temperature) using the methodology explained above. Take a printscreen of the scale, save it, and on the Image Color Picker webpage, choose the option “Use your Image”, tab “Upload Image” and upload your printscreen. Check the hex values and change the code. Do not forget to modify the maximum and minimum contour limits, according to the values of the scale used on the CPTEC webpage.

This is the expected result:

Depending on the zoom and the resolution of the model, sometimes the contours are not as smooth as we wanted. Let’s see how to smooth the contours.

SCRIPT 8: SMOOTHING THE CONTOURS

Let’s create our eight script. Create and save a file called “script_08.py”.

So far, the only manipulation we’ve done with our array is the conversion from K to °C. In this example script, we are going to smooth the contours of our plot with the scipy.ndimage.zoom function, which increases the size of our array “x” times, interpolating the data.

First, let’s learn how to check the dimensions of our data. We can do this with the shape command:

#-----------------------------------------------------------------------------------------------------------
# Convert to Celsius
tmtmp = tmtmp - 273.15

# Print the array dimensions
print("Array dimensions :", tmtmp.shape)
#-----------------------------------------------------------------------------------------------------------

Copy the contents of the previous script, adding the above command, run it and see the output in the terminal:

In the example above, we see that our data matrix for the selected region is 81 x 121 pixels.

Now, let’s reduce the plot region, to highlight its low resolution:

# Select the extent [min. lon, min. lat, max. lon, max. lat]
extent = [-85.0, 10.00, -60.00, 25.00]

Run the script again. See that the dimension of our array has decreased:

Seeing the plot of the reduced region, the low resolution of our data is evident, with “squared” contours:

Let’s apply the zoom command to smooth the contours, as follows:

#-----------------------------------------------------------------------------------------------------------
# Convert from K to °C
tmtmp = tmtmp - 273.15

print("\nArray dimensions before smoothing:")
print(tmtmp.shape)

# Smooth the contours
import scipy.ndimage
tmtmp = scipy.ndimage.zoom(tmtmp, 3)
lats = scipy.ndimage.zoom(lats, 3)
lons = scipy.ndimage.zoom(lons, 3)

print("Array dimensions after smoothing:")
print(tmtmp.shape)
#-----------------------------------------------------------------------------------------------------------

This should be our complete script at this point:

Run the “script_08.py” script with the following command:

python script_08.py

See that our array has tripled in size:

The following image will be displayed, with the smooth contours:

Now let’s see how to work with multiple files.

SCRIPT 9: WORKING WITH MULTIPLE FILES

In the previous scripts, we worked with a single file. In this example, we are going to process multiple files using a single script. Note: There are several ways to process multiple files, in this procedure we will see an example.

Access the NOMADS server and download files f003, f006, f009, f012, f015, f018, f021 and f024, following the procedures seen at the beginning of the pre-course activity, for the same day and run chosen for the f000 file downloaded in the first step. Don’t forget to choose the same options as before:

At this point, you should have 9 GRIB files in your working directory:

Let’s create our ninth script. Create and save a file called “script_09.py”.

First, let’s create via code the directory where our plots will be stored. The command below creates a directory if it doesn’t exist.

# Animation directory
dir = "Animation"; os.makedirs(dir, exist_ok=True)

To read several files in the same script, we will use a “for” loop, and in each cycle of the loop, we will read a different GRIB file. For this, we will create three variables, where we will define the start time, the end time and the loop interval.

# Data you want to process in the loop
hour_ini = 0   # Init time  
hour_end = 24  # End time
hour_int = 3   # Interval

In the example above, we indicate that we want to read the files from hour 0 to hour 24, with an interval of 3 hours.

Using the variables above, we create our “for” loop:

for hour in range(hour_ini, hour_end + 1, hour_int):

In the example above, we added 1 to the end hour, as the Python “for” loop disregards the last value in the loop (we don’t want to disregard the value “24”, so we added 1).

Already inside the “for” loop, we need to indicate which file to be processed in that cycle. We do this as follows:

#-----------------------------------------------------------------------------------------------------------
# INPE / CPTEC Training: NWP Data Processing With Python - Script 9: Working With Multiple Files
# 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 
import os                                  # Miscellaneous operating system interfaces
#----------------------------------------------------------------------------------------------------------- 

# Animation directory
dir = "Animation"; os.makedirs(dir, exist_ok=True)

# Select the extent [min. lon, min. lat, max. lon, max. lat]
extent = [-100.0, 0.00, -40.00, 40.00]

#-----------------------------------------------------------------------------------------------------------

# GRIB file name without the last three characters
path = ("gfs.t00z.pgrb2full.0p50.f")

# Data you want to process in the loop
hour_ini = 0   # Init time  
hour_end = 24  # End time
hour_int = 3   # Interval

#-----------------------------------------------------------------------------------------------------------

for hour in range(hour_ini, hour_end + 1, hour_int):

    # File to process
    grib = path + str(hour).zfill(3)
    
    # If the file exists
    if (os.path.exists(grib)):
        
        # Process the file
        print("\nProcessing file: ", grib)

        # Read the GRIB file
        grib = pygrib.open(grib)
        
        # Select the variable
        grb = grib.select(name='2 metre temperature')[0]

On line 23, we declare the variable “path”, with the name of the GRIB file, disregarding the last three characters (which represent the time).

In line 35, within the cycle, we concatenate the string “hour” with the string “path”, representing it with three characters (via the command “zfill”). On line 41, if the file exists, we proceed with the processing, in the same way we were doing before.

Inside the loop, also change the command that generates the images, as follows:

# Save the image
plt.savefig('Animation//image_loop_' + str(hour) + '.png', bbox_inches='tight', pad_inches=0, dpi=100)

This will save the images inside the “Animation” folder (created via code) and will remove excess whitespace from the plots (as we specified the image size as 10 x 6 inches and this is not necessarily the size of our data).

This should be our complete script at this point:

Run the “script_09.py” script with the following command:

python script_09.py

IMPORTANT: We are going to process 9 different files, so the processing time will be considerably longer. Follow the processing steps in the terminal.

After running the script, we will see the 9 images inside the “Animation” folder:

With this, let’s see our last example script of this pre-course activity, and generate an animation with the generated images.

SCRIPT 10: CREATING AN ANIMATION

In this example, we are going to create an animated gif using the imageio library.

import imageio        # Python interface to read and write a wide range of image data

The script is relatively simple compared to the others.

Let’s create our tenth script. Create and save a file called “script_10.py”.

Below is the complete script:

#-----------------------------------------------------------------------------------------------------------
# INPE / CPTEC Training: NWP Data Processing With Python - Script 10: Creating an Animation
# Author: Diego Souza
#-----------------------------------------------------------------------------------------------------------
import imageio.v2 as imageio       # Python interface to read and write a wide range of image data

# Images we want to include in the GIF
files = ['image_loop_0.png', 'image_loop_3.png', 'image_loop_6.png', 'image_loop_9.png', 'image_loop_12.png', 
         'image_loop_15.png', 'image_loop_18.png', 'image_loop_21.png', 'image_loop_24.png']

# Create the GIF
images = []
for file in files:
    images.append(imageio.imread('Animation//' + file))

# Save the GIF
imageio.mimsave('Animation//animation.gif', images, fps=1)

In line 5, we import the imageio library, in lines 8 and 9, we create a list with the images we want to add to the GIF (generated in the previous script). On lines 13 and 14, we add the images to the GIF. On line 17, we save the GIF, informing that we want 1 image per second in our animation.

Run the “script_10.py” script with the following command:

python script_10.py

When running it, a file called “animation.gif” will be created in the “Animation” folder:

Open and visualize the generated animation:

With this script we finish this pre-course activities! In this activity we looked at some of the basic scripts for manipulating NWP data. During the course, we will review these scripts and see new scripts with new operations!

ACCESS TO THE SCRIPTS OF THE PRE-COURSE ACTIVITY AND THE SCRIPTS SHOWN DURING THE VIRTUAL COURSE

Find the scripts seen in the pre-course activity at the link below:

https://github.com/diegormsouza/gnc-a-caribbean/tree/main/NWP

In the course we will see almost 30 scripts, which will be made available at the link above.

Last updated on March 20, 2023, 18:00 UTC

Please send your feedback answering the poll: