This IPython Notebook can be read as either walkthrough of some pygaarst examples, or as a tutorial on plotting satellite imagery, either as flat images or on a map. The plotting examples work perfectly fine without pygaarst as long as you have
If instead you have 1-dimensional LAT and LONG coordinate arays, the second point is typically achieved by making a meshgrid. Observe that latitudes correspond to the y-coordinate.
XLON_matrix, YLAT_matrix = numpy.meshgrid(LON_vector, LAT_vector)
This tutorial deals with satellite imagery and derived data, which is typically distributed in one of the following file formats:
These are far from being the only data formats widely used in satellite products. Microwave/SAR sensor processing software often produces proprietary binary data that first needs to be reformated with special tools.
pygaarst currently implements generic GeoTIFF support as well as HDF5 support specific to VIIRS data. The other formats listed above are in the works.
First, we import the required libraries.
%pylab inline
import os, os.path
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
Populating the interactive namespace from numpy and matplotlib
For modules from the pygaarst package, prefer "from pygaarst import [module]" syntax. Or "from pygaarst import [module] as [shortname]".
from pygaarst import raster
A VIIRS HDF5 file (name slightly shortened) is opened by passing a file path (absolute or relative) to the constructor of the HDF5 object provided by pygaarst.raster:
viirsh5 = 'tests/data/SVDNB_npp_d20130725_t2022157_b09030_cspp_dev.h5'
viirsband = raster.HDF5(viirsh5)
INFO:root:Opening tests/data/SVDNB_npp_d20130725_t2022157_b09030_cspp_dev.h5
The viirsband object now knows some of its properties:
print viirsband.bandname
print viirsband.dirname
VIIRS-DNB-SDR_All tests/data
The data are in viirsband.dataobj. We can pass in the path to the radiance data, using the bandname read from the file metadata. This returns a numpy array:
viirsdataset = viirsband.dataobj['All_Data/VIIRS-DNB-SDR_All/Radiance']
viirsdataset.shape
(768, 4064)
The viirsband object also knows the name of its georeference file. If it is placed in the same directory as the file, it can be read automatically from the metadata. Here it is:
print viirsband.dataobj.attrs['N_GEO_Ref'][0][0]
print viirsband.geofilepath
GDNBO_npp_d20130725_t2022157_e2023399_b09030_c20130725204748850631_cspp_dev.h5 tests/data/GDNBO_npp_d20130725_t2022157_e2023399_b09030_c20130725204748850631_cspp_dev.h5
viirsband.lats
array([[ 68.66014099, 68.66169739, 68.66325378, ..., 59.24265671, 59.23745346, 59.23223495], [ 68.66703033, 68.66858673, 68.67014313, ..., 59.24757385, 59.24237061, 59.2371521 ], [ 68.67391968, 68.67547607, 68.67702484, ..., 59.252491 , 59.24728775, 59.24206924], ..., [ 73.37992859, 73.38204193, 73.3841629 , ..., 62.3796196 , 62.37378693, 62.36794662], [ 73.38665009, 73.38877869, 73.39089203, ..., 62.38377762, 62.37794495, 62.37210464], [ 73.39338684, 73.39550781, 73.39762115, ..., 62.38793945, 62.38210678, 62.37627411]], dtype=float32)
Alternatively, we can provide the path to the georeference file when instantiating the viirsband object. It is expected to be in HDF5 format as well.
viirsband = raster.HDF5(viirsh5, geofilepath="[path/to/georeference/file.h5]")
The data object can be plotted as a simple image plot. However, the image will be flipped up/down and may be flipped left/right, depending on whether the path is ascending or descending. The image is 768 by 4064 pixels, a very wide strip.
fig = plt.figure(1, figsize=(18, 15))
plt.imshow(viirsdataset, cmap='bone', aspect='equal')
<matplotlib.image.AxesImage at 0x10831ddd0>
Alternatively, we can plot it in a lat/long coordinate system to see roughly where the data falls on the earth. As we can see, the image is not only warped (as expected), but also flipped both ways.
plt.figure(2, figsize=(18, 15))
plt.pcolormesh(viirsband.lons, viirsband.lats, viirsdataset, cmap='bone')
<matplotlib.collections.QuadMesh at 0x11abc5950>
For plotting on a map background, the data is subsetted. The bit that sticks into Alaska is at the end of the array, and we take the last 2000 columns. Then generate a matplotlib Basemap of Alaska (Albers Equal Area), transform the lat/long arrays to map coordinates, and plot the data on top with pyplot.pcolormesh().
viirsdata = viirsdataset[:, -2000:]
lat = viirsband.lats[:, -2000:]
lon = viirsband.lons[:, -2000:]
fig = plt.figure(3, figsize=(18, 15))
m = Basemap(width=1800000, height=1200000, resolution='i', projection='aea', lat_1=55., lat_2=75., lat_0=65., lon_0=-150.)
m.drawrivers(color='lightskyblue', zorder=20)
m.drawcoastlines(color='lightskyblue', zorder=20)
m.drawcoastlines(color='.08')
m.drawcountries(linewidth=1.5, color='chocolate')
m.fillcontinents(lake_color='lightskyblue', color='snow')
m.drawmapboundary(fill_color='lightskyblue')
# transform coordinates from lat/lon to map coordinates
xx, yy = m(lon, lat)
m.pcolormesh(xx, yy, viirsdata, cmap='bone', zorder=15)
m.drawmeridians(np.arange(-180, -50, 10), labels=[0, 0, 1, 1])
m.drawparallels(np.arange(30, 80, 5), labels=[1, 1, 0, 0])
{60: ([<matplotlib.lines.Line2D at 0x10b525c90>], [<matplotlib.text.Text at 0x10b528890>, <matplotlib.text.Text at 0x10b5289d0>]), 65: ([<matplotlib.lines.Line2D at 0x10b528090>], [<matplotlib.text.Text at 0x10b528910>, <matplotlib.text.Text at 0x10b528a90>]), 70: ([<matplotlib.lines.Line2D at 0x10b528450>], [])}
And a higher-resolution zoomed-in version:
fig = plt.figure(4, figsize=(18, 15))
m = Basemap(width=600000, height=400000, resolution='f', projection='aea', lat_1=55., lat_2=75., lat_0=65., lon_0=-150.)
m.drawrivers(color='lightskyblue', zorder=20)
m.drawcoastlines(color='lightskyblue', zorder=20)
m.drawcoastlines(color='.08')
m.drawcountries(linewidth=1.5, color='chocolate')
m.fillcontinents(lake_color='lightskyblue', color='snow')
m.drawmapboundary(fill_color='lightskyblue')
xx, yy = m(lon, lat)
m.pcolormesh(xx, yy, viirsdata, cmap='bone', zorder=15)
m.drawmeridians(np.arange(-180, -50, 5), labels=[0, 0, 1, 1])
m.drawparallels(np.arange(30, 80, 5), labels=[1, 1, 0, 0])
{65: ([<matplotlib.lines.Line2D at 0x10ad50e90>], [<matplotlib.text.Text at 0x10ad55310>, <matplotlib.text.Text at 0x10ad55390>])}
The pygaarst.raster module also provides a class to represent a generic GeoTIFF object. Below a few methods that can be applied to any GeoTIFF file. pygaarst.raster also provides a Landsatscene class (initialized with the directory path that contains all files of a scene) as well as a Landsatband class that inherits from GeoTIFF. The latter are aware of Landsat metadata and need to be documented elsewhere.
rgbgeotiff = os.path.join('tests/data/', 'LC80680152013194_754_crop2.tiff')
ls = raster.GeoTIFF(rgbgeotiff)
This file is a cropped subset of a Landsat 8 scene, converted to 8-bit for easier display, as a false color RGB using bands 7, 6, 4 (corresponding to 7, 5, 5 for previous landsat TM/ETM+ sensors). If we want to display the whole thing, the can be done with matplotlib after re-arranging :
print type(ls.data)
img = np.rollaxis(ls.data, 0, 3)
fig = plt.figure(5, figsize=(18, 15))
plt.imshow(img)
<type 'numpy.ndarray'>
<matplotlib.image.AxesImage at 0x10ad70bd0>
ls.proj4
'+proj=utm +zone=6 +datum=WGS84 +units=m +no_defs '
ls.projection
'PROJCS["WGS 84 / UTM zone 6N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-147],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32606"]]'
print ls.data.shape
print ls.filepath
(3, 793, 787) tests/data/LC80680152013194_754_crop2.tiff
ls.simpleplot() # plots each band in small
Show the first band (which happens to be a Landsat B7) alone. The bright spots are actively burning areas in the Stuart Creek 2 fire.
fig = plt.figure(6, figsize=(18, 15))
plt.imshow(ls.data[0,:,:], cmap='bone')
<matplotlib.image.AxesImage at 0x10ae87250>
If we use pcolormesh, the plot is squashed and upside-down. Careful.
fig = plt.figure(7, figsize=(18, 15))
plt.pcolormesh(ls.data[0,:,:], cmap='bone')
<matplotlib.collections.QuadMesh at 0x10abe1890>