VIIRS with Pytroll

The NOAA/NASA weather satellite Suomi National Polar Partnership (NPP) satellite was successfully launched October 28th, 2011.

Suomi NPP carries the Visible Infrared Imaging Radiometer Suite (VIIRS) which is a moderate resolution Imager with herritage from MODIS onboard the EOS satellites Aqua and Terra, the AVHRR onboard the NOAA and Metop platforms, and the OLS onboard the DMSP satellites. VIIRS has a nearly constant spatial resolution across a very large swath (width is around 3000km).

In preparation for VIIRS we made a plugin reader and custom compositer to mpop. This was tested on synthetic (sample) data prior to launch, allowing us to produce VIIRS imagery as soon as Direct Readout was turned on.

The example images below have been generated from direct readout data from the station in Norrköping. We have been using the Community Satellite Processing Package (CSPP) to go from RDR to SDR - see cspp

For this tutorial template config files (see Installation and configuration) can be used. These are located in the etc dir of the mpop source. Copy mpop.cfg.template, areas.def.template and npp.cfg.template to another dir and remove the .template extension. In the config file npp.cfg locate the section viirs-level2 and modify the defined dir to point to the directory where you keep the SDR data.

Set PPP_CONFIG_DIR to the directory containing your modified mpop config files.


>>> from mpop.satellites import PolarFactory
>>> from datetime import datetime
>>> time_slot = datetime(2013, 6, 25, 11, 15)
>>> orbit = "08599"
>>> global_data = PolarFactory.create_scene("npp", "", "viirs", time_slot, orbit)
>>> global_data.load([1.38])
>>> print global_data
 -------> print(global_data)
 'M01: (0.402,0.412,0.422)μm, resolution 742m, not loaded'
 'M02: (0.436,0.445,0.454)μm, resolution 742m, not loaded'
 'M03: (0.478,0.488,0.498)μm, resolution 742m, not loaded'
 'M04: (0.545,0.555,0.565)μm, resolution 742m, not loaded'
 'M05: (0.662,0.672,0.682)μm, resolution 742m, not loaded'
 'M06: (0.739,0.746,0.754)μm, resolution 742m, not loaded'
 'M07: (0.846,0.865,0.885)μm, resolution 742m, not loaded'
 'M08: (1.230,1.240,1.250)μm, resolution 742m, not loaded'
 'M09: (1.371,1.378,1.386)μm, shape (768, 3200), resolution 742m'
 'M10: (1.580,1.610,1.640)μm, resolution 742m, not loaded'
 'M11: (2.225,2.250,2.275)μm, resolution 742m, not loaded'
 'M12: (3.610,3.700,3.790)μm, resolution 742m, not loaded'
 'M13: (3.973,4.050,4.128)μm, resolution 742m, not loaded'
 'M14: (8.400,8.550,8.700)μm, resolution 742m, not loaded'
 'M15: (10.263,10.763,11.263)μm, resolution 742m, not loaded'
 'M16: (11.538,12.013,12.489)μm, resolution 742m, not loaded'
 'I01: (0.600,0.640,0.680)μm, resolution 371m, not loaded'
 'I02: (0.845,0.865,0.884)μm, resolution 371m, not loaded'
 'I03: (1.580,1.610,1.640)μm, resolution 371m, not loaded'
 'I04: (3.580,3.740,3.900)μm, resolution 371m, not loaded'
 'I05: (10.500,11.450,12.300)μm, resolution 371m, not loaded'
 'DNB: (0.500,0.700,0.900)μm, resolution 742m, not loaded'

We have now loaded the VIIRS M09 band. Now let us look at the data:

>>> img = global_data.image.channel_image(1.38)
>>> img.enhance(stretch='histogram')

The black stripes are due to the so called bowtie deletion, which is handled onboard the satellite. The bowtie effect is a geometric feature of the VIIRS scan. Similar to the MODIS sensor individual VIIRS lines will overlap as one approach the edge of the swath. These overlapping samples/pixels have been removed onboard in order to minimise the bandwidth usage on the broadcast. Thus when data are displayed un-projected these no-data lines will appear in the image.

Making RGB’s

Here is an example making a true color RGB with the VIIRS bands:

>>> global_data.load(global_data.image.truecolor.prerequisites)
>>> img = global_data.image.truecolor()

Some words on geolocation

Before we continue projecting the data to a map projection, just some words on geolocation. For both the M- and the I-bands there are normally two geolocation files available. The ones using a terrain correction are prefixed GMTCO for the M-bands and GITCO for the I-bands. Correspondingly there are two geolocation files where the topography has been neglected, and these are prefixed GMODO and GIMGO. For the Day-Night band there shall be only one geolocation file, prefixed GDNBO. In the configuration file you have some flexibility to constrain what type of geolocation to apply (terrain corrected or the non-corrected geoid files). However, normally it would be sufficient to allow all kinds of geolocations and let the viirs reader in mpop choose the appropriate one for you:

geo_filenames = G????_%(satellite)s_d%Y%m%d_t%H%M???_e???????_b%(orbit)s_c*_cspp_dev.h5
filename = SV???_%(satellite)s_d%Y%m%d_t%H%M???_e???????_b%(orbit)s_c*_cspp_dev.h5
dir = /data/viirs/sdr
format = viirs_sdr.ViirsSDRReader

If the list of geo-files is not specified in the config file or some of them do not exist, the viirs reader will rely on what is written in the band-data metadata header. For each band the SDR files contain an attribute giving the name of the matching geolocation file. The reader will assume that this gile is lying in the same directory as the band data given by the path dir above. If this file is not there no geolocation will be loaded.

With the above configuration, the viirs reader will search first for the terrain corrected files matching the string above. If the terrain corrected files are not available, the reader will try the geoid ones. If no geoid files matching the string above are found the reader will not load any geolocation. Without geolocation the data cannot be projected to a map-projetion, of course.

Map projection

Reprojecting data is done in exactly the same way the AVHRR data was reprojected in the Quickstart with AVHRR tutorial:

>>> local_data = global_data.project("scan500m", mode="nearest")
>>> img = local_data.image.truecolor()

Here we have defined an area called area500m covering Scandinavia, and with a pixel resolution of 500 meters. This definition is stored in the areas.def.template file. See the Quickstart with AVHRR tutorial.

It is easier to navigate in the image if we add coastlines and poltical boarders, so lets do that with PIL and pycoast:

>>> from PIL import Image
>>> from pycoast import ContourWriter
>>> from mpop.projector import get_area_def
>>> cw = ContourWriter('/local_disk/data/shapes')
>>> img ='./viirs_truecolor_proj.png')
>>> area_def = get_area_def("scan500m")
>>> cw.add_coastlines(img, area_def, resolution='i', level=3)

But what we actually wanted was to load all the available data (VIIRS granules) received that covers the area. To do this in a smart and economical way we only load the granules that are inside the area of interest. The VIIRS reader in mpop version 1.x and later supports loading of granules within a certain time interval. Optimally we would like to load granules by area coverage, but this is not yet possible in mpop alone. With the functionality provided by pyorbital and pyresample it will however be possible to do this. We leave this for an other time, and instead we try load the granules by specifying a time interval that we think will cover the area of interest:

>>> from mpop.satellites import PolarFactory
>>> from datetime import datetime
>>> orbit = "08599"
>>> tslot = datetime(2013, 6, 25, 11, 11)
>>> starttime = datetime(2013, 6, 25, 11, 11)
>>> endtime = datetime(2013, 6, 25, 11, 18)
>>> glbd = PolarFactory.create_scene("npp", "",
...                                  "viirs", tslot,
...                                   orbit)
>>> glbd.load(glbd.image.green_snow.prerequisites |
...           glbd.image.natural.prerequisites,
...           time_interval=(starttime, endtime))
>>> img = glbd.image.natural()

And now lets project it to the area:

>>> local_data = glb_data.project(areaid, mode="nearest", radius=2000)

We can display the green_snow composite as we already made sure to load the necessary channels earlier (see code above):

>>> img = local_data.image.green_snow()

High resolution images

The VIIRS sensor have 5 AVHRR-like channels with a resolution 3 times higher or even better (at edge of swath). These are the I-bands seen in the list above. Making imagery from these goes exactly the same way as for the M-bands. However, since there is overlap in the spectral range between I-bands and M-bands, you need to specify also the resolution or use the band name when loading:

>>> global_data.load(['I03'])
>>> global_data['I03'].show()

Generating and mapping the overview of the I-bands is done in the same way as for the M-bands of course. Here we have made a specific I-band overview method called hr_overview:

>>> from mpop.satellites import PolarFactory
>>> from datetime import datetime
>>> orbit = "08599"
>>> time_slot = datetime(2013, 6, 25, 11, 15)
>>> global_data = PolarFactory.create_scene("npp", "", "viirs", time_slot, orbit)
>>> global_data.load(global_data.image.hr_overview.prerequisites)
>>> local_data = global_data.project("scan500m", mode="nearest")
>>> img = local_data.image.hr_overview()

The Day/Night Band

The VIIRS Day/Night band draws heritage from the DMSP Operational Linescan System (OLS) and is a broad band channel in the Visible and Near-Infrared spectral range. It operates with three different gains to optimise the sensitivity independant of illumination. We find a nighttime case with some moonlight, and make a stretched black and white image for display:

>>> time_slot = datetime(2012, 8, 31, 1, 8)
>>> orbit = "04365"
>>> global_data = PolarFactory.create_scene("npp", "", "viirs", time_slot, orbit)
>>> global_data.load(['DNB'])
>>> global_data.image.dnb().show()

During nighttime it is sufficiently sensitive so that useful information on clouds and surfaces may be deduced from reflected moonlight. Naturally the units of this band cannot be given as a solar reflectance factor, but instead the radiance is provided:

>>> print global_data['DNB'].info
    {'units': 'W sr-1 m-2', 'band_id': 'DNB'}

The units in the HDF5 SDR file is W cm^{-2} sr^{-1} (see table 2.18.2-1, page 355 of the NPOESS Common Data Format Control Book - Volume III - D34862-03 Rev E CDRL No. A014). But in pytroll we keep to the physical units dictated by the netCDF CF convention on metadata, which is W m^{-2} sr^{-1}.

Observe that this is really the spectral radiance integrated over the entire band of wavelengths from 500 to 900 nm, and not a spectral radiance (e.g. unit W/(sr*m²*μm) which is otherwise common for narrow band channels.

>>> print global_data['DNB'].data
[[-- 0.000163395772688 0.000175373905222 ..., 0.000117216048238
  0.000115083799756 0.000105939005152]
 [-- 0.000170526851434 0.000165672157891 ..., 0.000114827875223
  0.000116445145977 0.000109282387712]
 [-- 0.000161146308528 0.000150438645505 ..., 0.000113979120215
  0.000117577423225 0.000114215945359]
 [-- 0.0001579762029 0.000168871047208 ..., 5.65401896893e-05
  5.81711428822e-05 6.36076947558e-05]
 [-- 0.000155943780555 0.000146388818393 ..., 6.1892089434e-05
  6.02728396188e-05 5.55949372938e-05]
 [-- 0.000157967660925 0.000157781381859 ..., 6.59923025523e-05
  6.20885184617e-05 5.79988991376e-05]]

We can check the range of radiaces in the granule and in print it in the units given in the input file if we like:

>>> print (global_data['DNB'].data * 10000).min()
>>> print (global_data['DNB'].data * 10000).max()

Let us load a few granules and assemble them and reproject them to get an image covering Scandinavia:

>>> time_slot = datetime(2013, 10, 20, 1, 26)
>>> starttime = datetime(2013, 10, 20, 1, 26)
>>> endtime = datetime(2013, 10, 20, 1, 34)
>>> orbit = "10252"
>>> global_data = PolarFactory.create_scene("npp", "", "viirs", time_slot, orbit)
>>> global_data.load(['DNB'], time_interval=(starttime, endtime))
>>> local_data = global_data.project('scan500m', mode="nearest", radius=2000)
>>> img = local_data.image.dnb(stretch='linear')