Blended TPW animation made with Python. The data is derived from several data sources, including the Special Sensor Microwave / Imager (SSM/I) passes from the Defense Meteorological Satellites Program (DMSP), the Advanced Microwave Sounding Unit (AMSU) from Polar Operational Environmental Satellites (POES), and information from Global Positional Satellite (GPS) meteorological data.
TPW Anomaly animation made with Python. The TPW anomaly assists forecasters to answer the question, “How unusual is the moisture field?”. Very high percentage values (200% or more) indicate a strong flooding potential or a possible severe weather indicator, while low values indicate potential fire hazards.
Hi GEONETCasters,
As seen on this and this Blog posts, you may find hourly Blended Total Precipitable Water and Total Precipitable Water Anomaly Products inside the JPSS/PRODUCTS/BTPW GNC-A Folder.
The file naming convention for these files are:
- NPR.COMP.TPW.Syydddhhmm.Eyydddhhmm.Thhmmss.he4
- NPR.COMP.PCT.Syydddhhmm.Eyydddhhmm.Thhmmss.he4
Where:
yy = year
jjj = julian day
hh = hour
mm = minutes
ss = seconds
D = day
S = start
E = end
T = time
These represent the year and julian day of the data covered, the start (“S”) and end (“E”) time in hour and minutes for which the data of the image covers, and the creation time (“T”) of the image in julian day, hour, and minutes.
Please find below an example Python script to manipulate these products.
Within the script you may also choose the region you want to plot. In the examples below the Brazilian region was chosen (BTPW on the left and TPW anomaly on the right).
In order to run the script, you need to install the pyhdf library. Here’s the conda command to download it:
conda install -c conda-forge pyhdf
Also, you must have the following text files referenced in the script (click to download):
###############################################################################
# LICENSE
#Copyright (C) 2018 - INPE - NATIONAL INSTITUTE FOR SPACE RESEARCH
#This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
#This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.
###############################################################################
# Required libraries ==========================================================
import matplotlib.pyplot as plt # Import the Matplotlib package
from mpl_toolkits.basemap import Basemap # Import the Basemap toolkit
import numpy as np # Import the Numpy package
import datetime # Library to convert julian day to dd-mm-yyyy
# conda install -c conda-forge pyhdf
from pyhdf.SD import SD, SDC
# =============================================================================
# Reading the file:
path = 'E:\\VLAB\\Python\\TPW Samples\\NPR.COMP.TPW.S180310337.E180311226.T134745.he4'
file = SD(path, SDC.READ)
#print (file.info())
datasets_dic = file.datasets()
#for idx,sds in enumerate(datasets_dic.keys()):
#print (idx,sds)
# Detecting if is BTPW or PCT:
product = (path[path.find("COMP.")+5:path.find(".S")])
if (product == "TPW"):
sds_obj = file.select('TPW_Mean') # select sds
product = "BTPW"
elif (product == "PCT"):
sds_obj = file.select('TPW_Percent') # select sds
product = "BTPW_PCT"
data = sds_obj.get() # get sds data
# The pixels without values are set with "-1". You may replace it with a nan if you want
#data[data == -1] = np.nan
# =============================================================================
# Reading Longitudes ==========================================================
with open('E:\\VLAB\\Python\\TPW Samples\\A_BTPW_Longitude.txt') as f:
longitudes = f.readlines()
num_linhas = len(longitudes)
lon = []
for x in range (2,num_linhas):
valores = longitudes[x].strip()
valores = valores.split()
lon.append(float(valores[2]))
f.close()
lons = np.asarray(lon)
# Reading Latitudes ===========================================================
with open('E:\\VLAB\\Python\\TPW Samples\\A_BTPW_Latitude.txt') as f:
latitudes = f.readlines()
num_linhas = len(latitudes)
lat = []
for x in range (2,num_linhas):
valores = latitudes[x].strip()
valores = valores.split()
lat.append(float(valores[2]))
f.close()
lats = np.asarray(lat)
# =============================================================================
# Defining the plot extent:
# GOES-East
#extent = [-156.29949951171875, -60.0, 6.299499988555908, 60.0]
# South America
#extent = [-90.0, -60.0, -30.0, 15.0]
# Full Extent
extent = [-180, -71.24783, 180, 71.24783]
# Latitude lower and upper index:
latli = np.argmin( np.abs( lats - extent[1] ) )
latui = np.argmin( np.abs( lats - extent[3] ) )
if (extent[2] > 20):
# Longitude lower and upper index:
lonli = np.argmin( np.abs( lons - extent[0] ) )
lonui = np.argmin( np.abs( extent[2] - lons ) )
data1 = data[latui:latli , lonli:2499]
data2 = data[latui:latli , 0:lonui]
data = np.hstack((data1, data2))
else:
# Longitude lower and upper index:
lonli = np.argmin( np.abs( lons - extent[0] ) )
lonui = np.argmin( np.abs( lons - extent[2] ) )
data = data[latui:latli , lonli:lonui]
# Define the size of the saved picture=========================================
DPI = 150
fig = plt.figure(figsize=(data.shape[1]/float(DPI), data.shape[0]/float(DPI)), frameon=False, dpi=DPI)
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)
ax = plt.axis('off')
#==============================================================================
# Create the Basemap
bmap = Basemap(projection='merc',llcrnrlon=extent[0], llcrnrlat=extent[1],urcrnrlon=extent[2], urcrnrlat=extent[3], resolution='l')
# Draw parallels and meridians
bmap.drawparallels(np.arange(-90.0, 90.0, 5.0), labels=[False,False,False,False], dashes=[5, 5], linewidth=0.2, fontsize=8, color='grey')
bmap.drawmeridians(np.arange(0.0, 360.0, 5.0), labels=[False,False,False,False], dashes=[5, 5], linewidth=0.2, fontsize=8, color='grey')
# Add the brazilian states shapefile
bmap.readshapefile('E:\\VLAB\\Python\\Shapefiles\\estados_2010','estados_2010',linewidth=0.30,color='gray')
# Add the countries shapefile
bmap.readshapefile('E:\\VLAB\Python\\Shapefiles\\ne_10m_admin_0_countries','ne_10m_admin_0_countries',linewidth=0.50,color='black')
# Add the continents shapefile
bmap.readshapefile('E:\\VLAB\\Python\\Shapefiles\\continent','continent',linewidth=1.00,color='black')
# Show image
if (product == "BTPW"):
bmap.imshow(data, origin='upper', cmap='rainbow', vmin = 0, vmax = 70)
elif (product == "BTPW_PCT"):
bmap.imshow(data, origin='upper', cmap='rainbow', vmin = 0, vmax = 300)
# Getting information from the file name ======================================
# Search for the Scan start in the file name
Start = (path[path.find(".S")+2:path.find(".E")])
# Converting from julian day to dd-mm-yyyy
year = int('20' + (Start[0:2]))
dayjulian = int(Start[2:5]) - 1 # Subtract 1 because the year starts at "0"
dayconventional = datetime.datetime(year,1,1) + datetime.timedelta(dayjulian) # Convert from julian to conventional
#print (dayconventional)
date = dayconventional.strftime('%d%m%Y') # Format the date according to the strftime directives
date_formatted = dayconventional.strftime('%d-%m-%Y')
# Search for time in the filename
Time = path[-10:]
Hour = Time [0:2]
Minutes = Time [2:4]
Seconds = Time [4:6]
time = Hour + ":" + Minutes + ":" + Seconds + " UTC"
# Colorbar ====================================================================
# Plot a colorbar
if (product == "BTPW"):
cb = bmap.colorbar(location='bottom', size = '2%', ticks=[5,10,15,20,25,30,35,40,45,50,55,60,65], pad = '-2.0%')
elif (product == "BTPW_PCT"):
cb = bmap.colorbar(location='bottom', size = '2%', ticks=[25,50,75,100,125,150,175,200,225,250,275], pad = '-2.0%')
# Remove the colorbar outline
cb.outline.set_visible(False)
# Remove the colorbar ticks
cb.ax.tick_params(width = 0)
# Put the colobar labels inside the colorbar
cb.ax.xaxis.set_tick_params(pad = -14)
# Change the color and size of the colorbar labels
cb.ax.tick_params(axis='x', colors='black', labelsize=10)
# =============================================================================
# Legend ======================================================================
if (product == "BTPW"):
cb.set_label('GNC-A Blog: Blended Total Precipitable Water [mm] ' + date_formatted + " " + time, color='white', size = 15)
elif (product == "BTPW_PCT"):
cb.set_label('GNC-A Blog: Total Precipitable Water Anomaly [%] ' + date_formatted + " " + time, color='white', size = 15)
cb.ax.xaxis.set_label_position('top')
cb.ax.xaxis.set_label_coords(0.245, 1.5)
# =============================================================================
# Add logos / images to the plot ==============================================
logo_GNC = plt.imread('E:\\VLAB\\Python\\Logos\\GNC Logo - White.png')
logo_INPE = plt.imread('E:\\VLAB\\Python\\Logos\\INPE Logo.png')
logo_NOAA = plt.imread('E:\\VLAB\\Python\\Logos\\NOAA Logo.png')
logo_JPSS = plt.imread('E:\\VLAB\\Python\\Logos\\JPSS Logo.png')
plt.figimage(logo_GNC, 10, 80, zorder=3, alpha = 1, origin = 'upper')
plt.figimage(logo_INPE, 390, 80, zorder=3, alpha = 1, origin = 'upper')
plt.figimage(logo_NOAA, 490, 80, zorder=3, alpha = 1, origin = 'upper')
plt.figimage(logo_JPSS, 590, 80, zorder=3, alpha = 1, origin = 'upper')
# =============================================================================
# =============================================================================
# Save the result
plt.savefig('E:\\VLAB\\Python\\Output\\JPSS_' + product + '_' + date + Hour + Minutes + Seconds + '.png', dpi=DPI, pad_inches=0)
plt.close()
For more details about the above files and what they contain, please visit the OSPO BTPW homepage http://www.ospo.noaa.gov/Products/bTPW/Overview.html, where you can find various kinds of information about the bTPW product and what it represents.
Please find below other blog posts about JPSS and GNC-A:
The top 10 visitor countries are:
GNC-A Blog 