This is the fourth part of the GOES-16 / Python tutorial series. You may find the other parts by clicking at the links below:

• Part 1: Extracting Values and Basic Plotting
• Part 2: Basemaps, Grids, Legend and Exporting the Result
• Part 3: Getting information from the file name and header

Today, we’ll learn how to add a background to our map and change the transparency of a layer.

Many satellite imagery viewers on the web, like Real Earth, overlays multiple images over a background map, giving the user the option to change the transparency and other features.

Putting our GOES-16 simulated channel over a background and changing the transparency is very simple with Python. By default, we have the following options for the background:

• Blue Marble
• Shaded Relief
• Etopo
• Land-Sea Mask

Let’s begin adding a Blue Marble background.

In order to do this, add the “bmap.bluemarble()” command before the “bmap.imshow” command. Also, on the “bmap.imshow” command, add the “alpha=0.5” parameter. The code below show the changes:

# Add a Blue Marble background
bmap.bluemarble()
# Plot GOES-16 channel with transparency
bmap.imshow(data, origin='upper', vmin=170, vmax=378, cmap='Greys', alpha=0.5)

Also, remove the legend from the code and add “Blue Marble Background” to the title.

By doing this, you should get this result:

We have our GOES-16 image with a 50% transparecy over a Blue Marble background. Isn’t this easy? Try changing the alpha parameter with different values in order to test different results.

Now try the same with the bmap.shadedrelief() and bmap.etopo() instructions. You should get these results:

Note: The Etopo plot was done with a 0.7 alpha.

Finally, we can also create ocean and land masks. To create a sea mask (very useful for land only products (like Land Surface Temperature), use the following command:

# Add an ocean mask
bmap.drawlsmask(ocean_color='aqua',lakes=True)

This should be the result:

If may choose any color for the ocean. At this page you may find some custom colors.

To add a land mask (very useful for products that cover only the ocean, like Sea Surface Temperature), you may use the following sintax:

# Add a land mask
bmap.drawlsmask(land_color='chartreuse')

This should be the result:

And this is the end of Part IV. Stay tuned for the next part!

This is the third part of the GOES-16 / Python tutorial series. You may find the other parts by clicking at the links below:

• Part 1: Extracting Values and Basic Plotting
• Part 2: Basemaps, Grids, Legend and Exporting the Result

Today, we’ll learn two things:

1-) How to extract the GOES-16 Band Number, the Start and the End of the Scan times based on the file name.

2-) How to read the GOES-16 NetCDF file header in order to retrieve information about the image.

By doing this, even if we change the source image, which in our example is declared at the following line …

# Path to the GOES-R simulated image file
path = 'C:\VLAB\OR_ABI-L2-CMIPF-M4C13_G16_s20161811455312_e20161811500135_c20161811500199.nc'

… the plot will automatically show the information of the new image. You’ll understand that when we finish this part.

Let’s start!

Extracting the Band Number, the Start and the End of the Scan times from the File Name:

In this blog post, we have seen that that the file naming convention for a GOES-16 channel image that will be on GNC-A follow this convention:

OR_ABI-L2–CMIPF–M4C13_G16_sYYYYJJJHHMMSSs_eYYYYJJJHHMMSSs_cYYYYJJJHHMMSSs.nc

Where:

OR: Operational System Real-Time Data
ABI-L2: Advanced Baseline Imager Level 2+
CMIPF: Cloud and Moisture Image Product – Full Disk
M3 / M4: ABI Mode 3 or ABI Mode 4
C13: Channel Number (Band 13 in this example)
G16: GOES-16
sYYYYJJJHHMMSSs: Observation Start
eYYYYJJJHHMMSSs: Observation End
cYYYYJJJHHMMSSs: File Creation

Note: YYYY = Year, JJJ = Julian Day, HH = Hour, MM = Month, SS = Seconds, s = Miliseconds

For now we will work with the Channel Number, Observarion Start and Observation End values.

What does Observation Start and Observation End means? Take a look at the video below and you will understand this concept:

Basically, the file name contains the time for the start of the globe scan and the end of the globe scan.

In our script, this is how the source image is specified:

# Path to the GOES-R simulated image file
path = 'C:\VLAB\OR_ABI-L2-CMIPF-M4C13_G16_s20161811455312_e20161811500135_c20161811500199.nc'

In order to extract the Band number, Observation Start and End, we need to search these values in the file name. In our case the file name is in the “path” variable. There are multiple ways to do this substring “search” with Python. This example will show you one way. Please add the following lines to the script, after creating the path variable:

# Search for the GOES-R channel in the file name
Band = (path[path.find("M4C")+3:path.find("_G16")])
print ("ABI Band: " + Band)

# Search for the Scan Start in the file name
Start = (path[path.find("s")+1:path.find("_e")])
print ("Observation Start: " + Start)

# Search for the Scan End in the file name
End = (path[path.find("e")+1:path.find("_c")])
print ("Observation End: " + End)

With this code we’re saying that we want three strings:

• The one between “M4C” and “_G16”, that is the channel
• The one between “s” and “_e”, that is the observation start
• The one between “e” and “_c”, that is the observation end.

When you execute the code, you should see the following output in the Spyder console:

The Python Console is located on the bottom right of the defult Spyder interface. You may choose the “IPython Console” tab (the image will be plotted inside the console) or the “Python Console” tab (a new window will be opened to show the image).

Note: The output location of the “print” command will depend on the Python development interface you are using.

Now add this line of code in the end of the script, after the “bmap.colorbar…” instruction and before the “DPI = 300…” instruction:

# Add a title to the plot
plt.title("GOES-16 ABI Simulated Band " + Band + " - Scan from " + Start + " to " + End)

This should be the result:

However, showing the times and dates as YYYYJJJHHMMSSs is not intuitive at all. Let’s arrange the times and dates in an easier way, adding these script lines after the creation of the Start and End variables:

# Format the "Observation Start" string
Start_Formatted = Start[0:4] + " Day " + Start[4:7] + " - " + Start [7:9] + ":" + Start [9:11] + ":" + Start [11:13] + "." + Start [13:14] + " UTC"
print (Start_Formatted)

# Format the "Observation End" string
End_Formatted = End[0:4] + " Day " + End[4:7] + " - " + End [7:9] + ":" + End [9:11] + ":" + End [11:13] + "." + End [13:14] + " UTC"
print (End_Formatted)

We’re simply separating parts of the “Start” string. When we use Start[0:4], we’re saying that we want the first for characters from the 20161811455312 string (the year).

When we use Start[4:7], we’re saying we want characters 5, 6 and 7 (the julian day), and so on.

Then, change the title line code to this:

# Add a title to the plot
plt.title("GOES-16 ABI Simulated Band " + Band + " - Scan from " + Start_Formatted + " to " + End_Formatted)

This should be the result:

Too much text in a single line? Just add a “\n” to the code and you will have multiple lines:

# Add a title to the plot
plt.title("GOES-16 ABI Simulated Band " + Band + "\n Scan from " + Start_Formatted + " to " + End_Formatted)

You should get this:

That’s much better! Now let’s see how to read the header from a NetCDF file.

Reading a NetCDF header with Python:

In order to know the available NetCDF header variables, just add the following lines to our previous code:

# Store the NetCDF file key variables in the "variable" variable, hahaha
variables = nc.variables.keys()
# Print the available variables
print (variables)

When you execute the code, you should see the following key variable list in the Spyder console:

As you may see, we have 36 key variables! Let’s read one of those, the “geospatial_lat_lon_extent”.

To read an especific key variable, use the following syntax:

# Read the header to retrieve the geospatial extent
geo_extent = nc.variables['geospatial_lat_lon_extent']
print (geo_extent)

This should be the output in the console:

You may see that inside the ‘geospatial_lat_lon_extent’ NetCDF variable, we have 11 sub variables (from “long_name” to “geospatial_lon_units”).

Reading an specific sub variable is very easy. The example line below shows how to read the center longitude:

geo_extent = nc.variables['geospatial_lat_lon_extent']
center = geo_extent.geospatial_lon_center print (center)

The output should be -89.5.

Following, this methodology, we can extract all the variables we want from a NetCDF file header, as below:

geo_extent = nc.variables['geospatial_lat_lon_extent']

# Extract the image bounds and center, converting to string
center = str(geo_extent.geospatial_lon_center)
west = str(geo_extent.geospatial_westbound_longitude)
east = str(geo_extent.geospatial_eastbound_longitude)
north = str(geo_extent.geospatial_northbound_latitude)
south = str(geo_extent.geospatial_southbound_latitude)

And use this information the way we want. Please note thar we converted the values to string, using str(), for we want to use these numbers as strings on the image.

In the example below we’ll put this information on the plot:

plt.text(-300000,300000,'Geospatial Extent \n' + west + '°W \n' + east + '°E \n' + north + '°N \n' + south + '°S \n' + 'Center = ' + center + '°', fontsize = 7)

This should be the result:

And this is the end of Part III!

For you’re reference, this is the complete script we have used in this part:

# GNC-A Blog Python Tutorial: Part III

# Required libraries
import matplotlib.pyplot as plt # Import the Matplotlib package
from netCDF4 import Dataset # Import the NetCDF Python interface
from mpl_toolkits.basemap import Basemap # Import the Basemap toolkit
import numpy as np # Import the Numpy package

# Path to the GOES-R simulated image file
path = 'C:\VLAB\OR_ABI-L2-CMIPF-M4C13_G16_s20161811455312_e20161811500135_c20161811500199.nc'
# Search for the GOES-R channel in the file name
Band = (path[path.find("M4C")+3:path.find("_G16")])
# Search for the Scan start in the file name
Start = (path[path.find("s")+1:path.find("_e")])
Start_Formatted = Start[0:4] + " Day " + Start[4:7] + " - " + Start [7:9] + ":" + Start [9:11] + ":" + Start [11:13] + "." + Start [13:14] + " UTC"
# Search for the Scan end in the file name
End = (path[path.find("e")+1:path.find("_c")])
End_Formatted = End[0:4] + " Day " + End[4:7] + " - " + End [7:9] + ":" + End [9:11] + ":" + End [11:13] + "." + End [13:14] + " UTC"

# Open the file using the NetCDF4 library
nc = Dataset(path)

# Extract the Brightness Temperature values from the NetCDF
data = nc.variables['CMI'][:]

# Create the basemap reference for the Satellite Projection
bmap = Basemap(projection='geos', lon_0=-89.5, lat_0=0.0, satellite_height=35786023.0, ellps='GRS80')

# Plot GOES-16 channel
bmap.imshow(data, origin='upper', vmin=170, vmax=378, cmap='Greys')

# Draw the coastlines, countries, parallels and meridians
bmap.drawcoastlines(linewidth=0.5, linestyle='solid', color='black')
bmap.drawcountries(linewidth=0.5, linestyle='solid', color='black')
bmap.drawparallels(np.arange(-90.0, 90.0, 10.0), linewidth=0.3, color='white')
bmap.drawmeridians(np.arange(0.0, 360.0, 10.0), linewidth=0.3, color='white')

# Insert the legend at the bottom
bmap.colorbar(location='bottom', label='Brightness Temperature [K]')

# Add a title to the plot
plt.title("GOES-16 ABI Simulated Band " + Band + "\n Scan from " + Start_Formatted + " to " + End_Formatted)

# Read some variables from the NetCDF header in order to use it in the plot
geo_extent = nc.variables['geospatial_lat_lon_extent']

center = str(geo_extent.geospatial_lon_center)
west = str(geo_extent.geospatial_westbound_longitude)
east = str(geo_extent.geospatial_eastbound_longitude)
north = str(geo_extent.geospatial_northbound_latitude)
south = str(geo_extent.geospatial_southbound_latitude)

# Put the information retrieved from the header in the final image
plt.text(-300000,300000,'Geospatial Extent \n' + west + '°W \n' + east + '°E \n' + north + '°N \n' + south + '°S \n' + 'Center = ' + center + '°', fontsize = 7)

# Export result
DPI = 300
plt.savefig('C:\VLAB\Channel_13_python.png', dpi=DPI, bbox_inches='tight', pad_inches=0)

# Show the plot
plt.show()

Stay tuned for the next tutorial!

In the first part of this tutorial series, we have learned how to extract the brightness temperature values and make a simple plot of a simulated GOES-16 channel using python and the NetCDF library.

In this part, we’re going to learn:

• How to apply basemaps, grids and a legend;
• How to save the plot in a PNG file.

These are the necessary steps:

1-) First of all, we need to install the Basemap Toolkit. In order to do this, open the Windows command prompt or the Linux terminal and insert and execute the following command:

conda install -c conda-forge basemap=1.1.0

Note: You need to choose Yes (‘y’) + Enter when asked if you want to proceed with the installation.

Update the Basemap Toolkit using the following command:

conda install -c conda-forge basemap-data-hires=1.0.8.dev0

2-) Paste the following code in the Spyder Editor:

# GNC-A Blog GOES-16 Python Tutorial: Part II

# Required libraries
import matplotlib.pyplot as plt # Import the Matplotlib package
from netCDF4 import Dataset # Import the NetCDF Python interface
from mpl_toolkits.basemap import Basemap # Import the Basemap toolkit
import numpy as np # Import the Numpy package

# Path to the GOES-R simulated image file
path = 'C:\VLAB\OR_ABI-L2-CMIPF-M4C13_G16_s20161811455312_e20161811500135_c20161811500199.nc'

# Open the file using the NetCDF4 library
nc = Dataset(path)

# Extract the Brightness Temperature values from the NetCDF
data = nc.variables['CMI'][:] 

# Create the basemap reference for the Satellite Projection
bmap = Basemap(projection='geos', lon_0=-89.5, lat_0=0.0, satellite_height=35786023.0, ellps='GRS80')

# Plot GOES-16 Channel using 170 and 378 as the temperature thresholds
bmap.imshow(data, origin='upper', vmin=170, vmax=378, cmap='Greys')

# Draw the coastlines, countries, parallels and meridians
bmap.drawcoastlines(linewidth=0.3, linestyle='solid', color='black')
bmap.drawcountries(linewidth=0.3, linestyle='solid', color='black')
bmap.drawparallels(np.arange(-90.0, 90.0, 10.0), linewidth=0.1, color='white')
bmap.drawmeridians(np.arange(0.0, 360.0, 10.0), linewidth=0.1, color='white')

# Insert the legend
bmap.colorbar(location='bottom', label='Brightness Temperature [K]')

# Export result
DPI = 300
plt.savefig('C:\VLAB\GOES-16_Ch13.png', dpi=DPI, bbox_inches='tight', pad_inches=0)

# Show the plot
plt.show()

UPDATE 1, on Sept. 28 2017: You may add the “resolution” parameter to the Basemap instruction, so the basemap resolution is set, as below:

bmap = Basemap(projection=’geos’, lon_0=-89.5, lat_0=0.0, satellite_height=35786023.0, ellps=’GRS80′, resolution = ‘i’).

The following options may be used: ‘c’ (crude), ‘l’ (low), ‘i’ (intermediate), ‘h’ (high) and ‘f’ (full). However, the higher the resolution, higher the processing time required.

UPDATE 2, on Sept. 28 2017: You may have get a high resolution PNG adding the following lines before the first bmap instruction:

DPI = 150
ax = plt.figure(figsize=(2000/float(DPI), 2000/float(DPI)), frameon=True, dpi=DPI)

Then , delete the “DPI = 300” in the end of the script.

I just noticed that on Part VII, so you may use this for now on.

3-) Execute the program by clicking at the “Play” icon in the upper menu or hitting “F5”. This should be the result:

And this is it!

Let’s explain the code we have added, in comparison to the first part.

With the following lines we have added the Basemap and Numpy libraries:

from mpl_toolkits.basemap import Basemap # Import the Basemap toolkit
import numpy as np # Import the Numpy package

With the following line we have created the Basemap reference with the satellite projection (we’re going to see other projections on the next tutorials).

bmap = Basemap(projection='geos', lon_0=-89.5, lat_0=0.0, satellite_height=35786023.0, ellps='GRS80')

For now, it’s important to understand that -89.5 is the current position of GOES-16. In the end of the year it will go to either East (-75) or West (-135) position. That’s the only number you’ll have to change in this line when that happens.

Then, with the following line we create the plot with a color gradient.

bmap.imshow(data, origin='upper', vmin=170, vmax=378, cmap='Greys')

Please not that the grayscale gradient is setted according to the minimum and maximum thresholds we have chosen (170 ~ 378 K). You may change this line with different thresholds and see the results.

If you remove the thresholds from the code, as this:

bmap.imshow(data, origin='upper', cmap='Greys')

The plot will automatically be adjusted to the minimum and maximum Brightness Temperature values from the image. This is not recommended because we need color consistency between different images (because we have different min’s and max’s for each image) .

In the cmap parameter, we have chosen the “Greys” option. You may choose any option from this link.

In the example below, the “jet” colormap and the 170 ~ 325 K thresholds are used.

The next lines configure the basemap, parallels and meridians appearance:

bmap.drawcoastlines(linewidth=0.3, linestyle='solid', color='black')
bmap.drawcountries(linewidth=0.3, linestyle='solid', color='black')
bmap.drawparallels(np.arange(-90.0, 90.0, 10.0), linewidth=0.1, color='white')
bmap.drawmeridians(np.arange(0.0, 360.0, 10.0), linewidth=0.1, color='white')

You may test different line widths and colors. Here’s an example with the “hsv” colormap, a threshold between 170 and 360, countries and coastlines with a line width of 0.5 and the white color, and the parallels and meridians with a line width of 0.3 and the black color:

With the next line of code we have chosen the legend position and text:

bmap.colorbar(location='bottom', label='Brightness Temperature [K]')

Here’s an example with the legend on the right position (just change to location=’right’):

Could you tell which colormap we used in this example? 🙂

Finally, with the last lines of code, we save the result in a PNG file called “GOES-16_Ch13.png” with 300 DPI of resolution:

DPI = 300
plt.savefig('C:\VLAB\GOES-16_Ch13.png', dpi=DPI, bbox_inches='tight', pad_inches=0)

All right! In the next tutorial we’ll see how to read the NetCDF file name automatically in order to get the time and dates and add this info to the image!

See you there!