GEONETClass: Manipulating GOES-16 Data With Python – Part XI

Python-NetCDF-11.png

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

In this part, we’ll learn:

  • How to generate georeferenced PNG’s (and other formats) in order to open the plots with Geographical Information Systems like SIGMACast (and others)  

-THE PROBLEM-

Until now, we have seen how to generate GeoTIFF’s (Part VII), and PNG’s. The PNG’s we have generated until now are not georeferenced. We cannot add these PNG’s as a layer on a Geographical Information System (GIS) with other data for example.

Let’s see how to do it, step by step. This will open our application range greatly.

-THE SOLUTION-

To georeference an image (PNG, JPEG, etc), apart to create the image itself, as we have been doing until now, we need to create an additional file, called the image “navigation” file.

The navigation file is very simple. It is just a plain text file with six lines, as follows:

  • Line 1: Image spacing in the X direction
  • Line 2: Y rotation
  • Line 3: X rotation
  • Line 4: Image spacing in the Y direction (negative)
  • Line 5: Upper-left corner coordinate in X.
  • Line 6: Upper-right corner coordinate in Y.

Apart from this, we have three requisites for the navigation file and image:

1-) The name of the navigation file must be the same of the image.

2-) The extension of the navigation file must follow a rule:

The extension must follow this pattern: first letter of the original extension + last letter of the original extension + the letter “w”.

For example, if we have an image with the “.jpg” extension, the navigation text file must have the “.jgw” extension. If we have an image with the “.png” extension, the navigation text file must have the “.pgw” extension. If we have an image with the “.tif” extension, the navigation text file must have the “.tfw” extension.

3-) The generated image must be in the equidistant projection, as we are doing since PART VI, thanks to the glorious GDAL 🙂

-A PRACTICAL EXAMPLE-

Let’s plot a GOES-16 image as we have been doing until now. Please download the NetCDF from this link. It is a GOES-16 NetCDF for Channel 13, March 2nd, 2018, 13:00 UTC.

Use the following extent, for South America:

# Choose the visualization extent (min lon, min lat, max lon, max lat)
extent = [-90.0, -60.0, -30.0, 15.0]

Let’s call the image G16_SA_CH13_201803021300.png (SA being “South America for example), adapting the following piece of code:

# Date as string
date_save = dayconventional.strftime('%d%m%Y')
# Time (UTC) as string
time_save = Start [7:9] + Start [9:11]
# Save the result
plt.savefig('E:\\VLAB\\Python\\Output\\G16_SA_C' + str(Band) + '_' + date_save + time_save + '.png', dpi=DPI, pad_inches=0, transparent=True)
plt.close()

This should be the result:

South_America_Plot.png

Now, let’s remove the shapefiles, parallels, meridians, the logos, the legend and the black rectangle with the image description. You probably not want those in a GIS, specially the shapefile.

After removing everything, your code should look like this:

###############################################################################
# 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/.
###############################################################################
#======================================================================================================
# GNC-A Blog Python Tutorial: Part XI
#======================================================================================================
# 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
from remap import remap # Import the Remap function
from cpt_convert import loadCPT # Import the CPT convert function
from matplotlib.colors import LinearSegmentedColormap # Linear interpolation for color maps
import datetime # Library to convert julian day to dd-mm-yyyy
from matplotlib.patches import Rectangle # Library to draw rectangles on the plot
from netCDF4 import Dataset # Import the NetCDF Python interface
import sys # Import the "system specific parameters and functions" module
#======================================================================================================
# Load the Data =======================================================================================
# Path to the GOES-16 image file
#path = sys.argv[1]
path = 'E:\\VLAB\\Python\\GOES-16 Samples - AWS\\OR_ABI-L2-CMIPF-M3C13_G16_s20180611300422_e20180611311200_c20180611315527.nc'
# Getting information from the file name ==============================================================
# Search for the Scan start in the file name
Start = (path[path.find("_s")+2:path.find("_e")])
# Search for the GOES-16 channel in the file name
Band = int((path[path.find("M3C" or "M4C")+3:path.find("_G16")]))
# Create a GOES-16 Bands string array
Wavelenghts = ['[]','[0.47 μm]','[0.64 μm]','[0.865 μm]','[1.378 μm]','[1.61 μm]','[2.25 μm]','[3.90 μm]','[6.19 μm]','[6.95 μm]','[7.34 μm]','[8.50 μm]','[9.61 μm]','[10.35 μm]','[11.20 μm]','[12.30 μm]','[13.30 μm]']

# Converting from julian day to dd-mm-yyyy
year = int(Start[0:4])
dayjulian = int(Start[4:7]) - 1 # Subtract 1 because the year starts at "0"
dayconventional = datetime.datetime(year,1,1) + datetime.timedelta(dayjulian) # Convert from julian to conventional
date = dayconventional.strftime('%d-%b-%Y') # Format the date according to the strftime directives

time = Start [7:9] + ":" + Start [9:11] + ":" + Start [11:13] + " UTC" # Time of the Start of the Scan

# Get the unit based on the channel. If channels 1 trough 6 is Albedo. If channels 7 to 16 is BT.
if Band <= 6:
Unit = "Reflectance"
else:
Unit = "Brightness Temperature [°C]"

# Choose a title for the plot
Title = " GOES-16 ABI CMI Band " + str(Band) + " " + Wavelenghts[int(Band)] + " " + Unit + " " + date + " " + time
# Insert the institution name
Institution = "GNC-A Blog"
# =====================================================================================================

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

# Get the latitude and longitude image bounds
geo_extent = nc.variables['geospatial_lat_lon_extent']
min_lon = float(geo_extent.geospatial_westbound_longitude)
max_lon = float(geo_extent.geospatial_eastbound_longitude)
min_lat = float(geo_extent.geospatial_southbound_latitude)
max_lat = float(geo_extent.geospatial_northbound_latitude)

# Choose the visualization extent (min lon, min lat, max lon, max lat)
# South America
extent = [-90.0, -60.0, -30.0, 15.0]
# Full Disk
#extent = [min_lon, min_lat, max_lon, max_lat]

# Choose the image resolution (the higher the number the faster the processing is)
resolution = 2

# Calculate the image extent required for the reprojection
H = nc.variables['goes_imager_projection'].perspective_point_height
x1 = nc.variables['x_image_bounds'][0] * H
x2 = nc.variables['x_image_bounds'][1] * H
y1 = nc.variables['y_image_bounds'][1] * H
y2 = nc.variables['y_image_bounds'][0] * H

# Call the reprojection funcion
grid = remap(path, extent, resolution, x1, y1, x2, y2)

# Read the data returned by the function
if Band <= 6:
data = grid.ReadAsArray()
else:
# If it is an IR channel subtract 273.15 to convert to ° Celsius
data = grid.ReadAsArray() - 273.15
# Make pixels outside the footprint invisible
data[data <= -180] = np.nan

#======================================================================================================

# Define the size of the saved picture=================================================================
#print (data.shape)
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')
#======================================================================================================

# Plot the Data =======================================================================================
# Create the basemap reference for the Rectangular Projection
bmap = Basemap(llcrnrlon=extent[0], llcrnrlat=extent[1], urcrnrlon=extent[2], urcrnrlat=extent[3], epsg=4326)

if Band  7 and Band  10:
# Converts a CPT file to be used in Python
cpt = loadCPT('E:\\VLAB\\Python\\Colortables\\IR4AVHRR6.cpt')
# Makes a linear interpolation
cpt_convert = LinearSegmentedColormap('cpt', cpt)
# Plot the GOES-16 channel with the converted CPT colors (you may alter the min and max to match your preference)
bmap.imshow(data, origin='upper', cmap=cpt_convert, vmin=-103, vmax=84)
# Insert the colorbar at the bottom

# Date as string
date_save = dayconventional.strftime('%d%m%Y')
# Time (UTC) as string
time_save = Start [7:9] + Start [9:11]

# Save the result
plt.savefig('E:\\VLAB\\Python\\Output\\G16_SA_C' + str(Band) + '_' + date_save + time_save + '.png', dpi=DPI, pad_inches=0)#, transparent=True)
plt.close()

with open('E:\\VLAB\\Python\\Output\\G16_Log.txt', 'a') as log:
log.write(path.replace('\\\\', '\\') + '\n')

Plot the image again. This should be the result:

South_America_Plot_No_Additions.png

You may open the PNG in any image viewer, or send this result to a webserver, like this one, from one of our readers. But you still can’t use it on a Geographical Information System (GIS) because you do not have the navigation file.

Let’s create it! Remember that we need a text file with the following lines:

  • Line 1: Image spacing in the X direction
  • Line 2: Y rotation (we don’t want to rotate, so we’ll keep it zero)
  • Line 3: X rotation (we don’t want to rotate, so we’ll keep it zero)
  • Line 4: Image spacing in the Y direction (negative)
  • Line 5: Upper-left corner coordinate in X.
  • Line 6: Upper-right corner coordinate in Y.

Add the following code in the end of your script:

# Create the georeference file

# Line 1: Image spacing in the X direction
resl_x = abs(extent[3]-extent[1]) / data.shape[0]
# Line 2 and 3 are zero! We don't want to rotate the image
# Line 4: Image spacing in the Y direction (negative)
resl_y = abs(extent[2]-extent[0]) / data.shape[1] * (-1)
# Line 5: Upper-left corner coordinate in X
lonO = extent[0]
# Line 6: Upper-right corner coordinate in Y
latN = extent[3]

# Format the text, one parameter per line
navigation = '{}\n{}\n{}\n{}\n{}\n{}'.format(resl_x, 0, 0, resl_y, lonO, latN)

# Write the file with the same file name as the image, and the ".pgw" extent
open('E:\\VLAB\\Python\\Output\\G16_SA_C' + str(Band) + '_' + date_save + time_save + '.pgw', 'w').write(navigation)

After executing the script, now you should see two files in your “output” folder, one with the “.png” extension, and the other with the “.pgw” extension, our navigation file.

South_America_Plot_With_PGW

Open the “G16_SA_C13_020320181300.pgw” file with a text editor. This is what you should see:

Pgw file content

This is exactly what we needed.

-BEHOLD THE MAGIC-

Open the PNG image file with a GIS that supports reading georeferenced imagery. Please find below how it seems when opened with QGIS.

South_America_Quantum.png

You may zoom it to any region, after all, you have the full resolution right? Of course, the higher the resolution, the higher your computer specs should be for the software to run smoothly. If you’re not able to see it smoothly, increase the “resolution” variable in the code.

South_America_Quantum_Zoom.png

Now, let’s see how it looks when opened with the upcoming SIGMACast software. After adding this product as a SIGMACast layer (we’ll see how to do it in the near future):

SIGMACast Example 8 Part XI

This is what it looks when opened.

SIGMACast Example 1 Part XI.png

You may also zoom it to any region:

SIGMACast Example 2 Part XI.png

… Add custom shapefiles:

SIGMACast Example 3 Part XI.png

… Increase the layer transparency:

SIGMACast Example 4 Part XI.png

… Reduce the transparency:

SIGMACast Example 5 Part XI.png

… Measure distance between points:

SIGMACast Distance.png

… See the image in the “geos” projection:

SIGMACast Example 6 Part XI.png

… Animate, make the interface available over the local network, among other things… But that’s for the future. We’ll keep you informed.

-CONCLUSION-

We have learned today how to create navigation files for our plots. As said in the beggining, this could open our application range greatly. We’ll study more about it when SIGMACast V3.0 is publicly released. Stay tuned!

Advertisements

One thought on “GEONETClass: Manipulating GOES-16 Data With Python – Part XI

  1. Hello,

    Thanks for your great website. I have a problem running the program.

    When I run the program, I get the following error with the remap function:

    AttributeError: ‘NoneType’ object has no attribute ‘SetProjection’ which is with the raw file (seems it doesn’t get any value). Can you please help me on that?

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s