#!/usr/bin/env python
# encoding: utf-8
"""
plot_viirs_RGB_crefl.py

The purpose of this script is to create an RGB image from several  
VIIRS bands, which have been atmospherically corrected by crefl.

###########
True-color
###########

>>>>>>>>>>>
Long Answer
>>>>>>>>>>>

True-color imagery uses MODIS Bands 1, 4, and 3, respectively corresponding 
to the red, green, and blue range of the light spectrum, are assigned to the 
red, green, and blue channels of a digital image. These images are called 
true-color or natural color because this combination of wavelengths is similar 
to what the human eye would see.
	
>>>>>>>>>>>>
Short Answer
>>>>>>>>>>>>

>>> Modis Bands: 1,4,3 (670 nm: 565 nm: 479 nm)
>>> Viirs Bands: M05,M04,M03 (672 nm, 555 nm ,488 nm)

>>> Colors: Fire detection=red outlines
>>> Advantages: Natural-looking images of land surface, oceanic and atmospheric features. 

######################
Band 3,6,7 Combination
######################

>>>>>>>>>>>
Long Answer
>>>>>>>>>>>

The 3-6-7 composite assigns Bands 3, 6, and 7 to the red, green, and blue components 
of a digital image. This combination is good for revealing snow and ice because they 
are very reflective in the visible part of the spectrum, and very absorbent in 
Bands 6 and 7, which are a part of the spectrum called the short-wave infrared, or SWIR.

>>> Snow and Ice
Since the only visible light used in these images (Band 3) is assigned to red, snow and 
ice appear bright red. The more ice, the stronger the absorption in the SWIR bands, 
and the more red the color. Thick ice and snow appear vivid red (or red-orange), while 
small ice crystals in high-level clouds will appear reddish-orange or peach.

>>> Vegetation
Vegetation is absorbent in Band 3 and Band 7, but reflective in Band 6, and so will appear 
greenish in this band combination. Bare soil will appear bright cyan in the image since it 
much more reflective in Band 6 and Band 7 than Band 3.

>>> Water
Liquid water on the ground will be very dark since it absorbs in the red and the SWIR, but 
small liquid water drops in clouds scatter light equally in both the visible and the SWIR, 
and will therefore appear white. Sediments in water appear dark red.
	
>>>>>>>>>>>>
Short Answer
>>>>>>>>>>>>

>>> Modis Bands: 3,6,7 (479 nm, 1652 nm, 2155 nm)
>>> Viirs Bands: M03, M10, M11 (488 nm, 1610 nm, 2250 nm)

>>> Colors: Vegetation=Green
        Ice or snow=Red
        Liquid water on the ground=Black or dark red
        Liquid water clouds=White
        Ice clouds=Peach
        Desert=Light blue-green

>>> Advantages: Distinguishing liquid water from frozen water, for example, clouds over snow, 
            ice cloud versus water cloud; or floods from dense vegetation. 

######################
Band 7-2-1 Combination
######################

>>>>>>>>>>>
Long Answer
>>>>>>>>>>>

In this composite, MODIS Bands 7, 2, and 1, are assigned to the red, green, and blue 
portions of the digital image. This combination is most useful for identifying burn scars.

>>> Vegetation and bare ground
Vegetation is very reflective in the near infrared (Band 2), and absorbent in Band 1 and 
Band 7. Assigning that band to green means even the smallest hint of vegetation will appear 
bright green in the image. Naturally bare soil, like a desert, is reflective in all bands 
used in this image, but more so in the SWIR (Band 7, red) and so soils will often have a 
pinkish tinge.

>>> Burned areas
If vegetation burns, bare soil in the area will become exposed. Band 1 slightly increases 
usually, but that may be offset by the presence of black carbon residue. The near infrared 
(Band 2) will become darker, and Band 7 becomes more reflective. When assigned to red in 
the image, Band 7 will show burn scars as deep or bright red, depending on the type of 
vegetation burned, the amount of residue, or the completeness of the burn.

>>> Water
As with the 3-6-7 composite, water will appear black. Sediments in water appear dark blue.
	
>>>>>>>>>>>>
Short Answer
>>>>>>>>>>>>

>>> Modis Bands: 7,2,1 (2155 nm, 876 nm, 670 nm)
>>> Viirs Bands: M11, M07, M05 (2250 nm, 865 nm, 672 nm)

>>> Colors: Vegetation=Green
        Water=Black or dark blue
        Desert/Naturally bare soil=Sandy pink
        Burn scar=Red to reddish-brown, depending on the nature of the pre-fire vegetation and the severity of the burn.

>>> Advantages
Distinguishing burn scars from naturally low vegetation or bare soil. Enhancing floods. 

Minimum commandline...

python plot_viirs_RGB_crefl.py --geo_file=geofile.h5 --red_file=redfile.h5 --green_file=greenfile.h5 --blue_file=bluefile.h5

     OR

python plot_viirs_RGB_crefl.py --geo_file=geofile.h5 --crefl_file=creflFile.h5

Created by Geoff Cureton on 2012-10-03.
Copyright (c) 2012 University of Wisconsin SSEC. All rights reserved.

"""

file_Date = '$Date$'
file_Revision = '$Revision$'
file_Author = '$Author$'
file_HeadURL = '$HeadURL$'
file_Id = '$Id$'

__author__ = 'G.P. Cureton <geoff.cureton@ssec.wisc.edu>'
__version__ = '$Id$'
__docformat__ = 'Epytext'

import os, sys, logging, traceback
from os import path,uname,environ
import string
import re
import uuid
import shlex, subprocess
from subprocess import CalledProcessError, call
from shutil import rmtree
from glob import glob
import optparse as optparse
from time import time
from datetime import datetime,timedelta

import numpy as np
from  numpy import ma as ma
import scipy as scipy

import matplotlib
import matplotlib.cm as cm
from matplotlib.colors import ListedColormap
from matplotlib.figure import Figure

matplotlib.use('Agg')
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas

#matplotlib.use('WXAgg')
#from matplotlib.backends.backend_wxagg import FigureCanvasAgg as FigureCanvas

# This must come *after* the backend is specified.
import matplotlib.pyplot as ppl

from mpl_toolkits.basemap import Basemap

import Image as Image

import VIIRS as VIIRS

import tables as pytables
from tables import exceptions as pyEx

from pyhdf.SD import SD,SDC,SDim
import pyhdf.HDF as pyHDF

# every module should have a LOG object
sourcename= file_Id.split(" ")
LOG = logging.getLogger(sourcename[1])

###################################################
#                  Global Data                    #
###################################################

#global someGlobalVar 

#def someFunction(someParameter) :
    #global someGlobalVar
    #someVar = someGlobalVar + someParemeter

###################################
#          Get File Lists         #
###################################

def granuleFiles(geoGlob,prodGlob):
    '''
    Returns sorted lists of the geolocation and product files.
    '''
    
    geoDir = path.dirname(path.abspath(path.expanduser(geoGlob)))
    prodDir = path.dirname(path.abspath(path.expanduser(prodGlob)))

    print "Initial geoGlob = ",geoGlob
    print "Initial prodGlob = ",prodGlob

    geoGlob = path.basename(path.abspath(path.expanduser(geoGlob)))
    prodGlob = path.basename(path.abspath(path.expanduser(prodGlob)))

    geoPrefix = string.split(geoGlob,'_')[0]

    print "geoDir = ",geoDir
    print "prodDir = ",prodDir
    print "geoGlob = ",geoGlob
    print "prodGlob = ",prodGlob
    print "geoPrefix = ",geoPrefix

    geoList_in = glob("%s/%s" % (geoDir,geoGlob))
    prodList_in = glob("%s/%s" % (prodDir,prodGlob))
    geoList_in.sort()
    prodList_in.sort()
    
    #print "prodList_in..."
    #for prodFile in prodList_in:
        #print prodFile


    geoList = []
    prodList = []
    #prodList = prodList_in
    for files in prodList_in :
        prod_arr = string.split(path.basename(files),"_")
        #print "prod_arr = ",prod_arr
        dateStamp = prod_arr[2]
        timeStamp = prod_arr[3]
        geoFileGlob="%s/%s*%s_%s*.h5" % (geoDir,geoPrefix,dateStamp,timeStamp)
        #print "dateStamp = ",dateStamp
        #print "timeStamp = ",timeStamp
        #print "geoFileGlob = ",geoFileGlob
        geoFile = glob("%s/%s*%s_%s*.h5" % (geoDir,geoPrefix,dateStamp,timeStamp))
        if (np.shape(geoFile)[0] != 0) :
            geoList.append(geoFile[0])
            prodList.append(files)
        else :
            #geoList.append(files)
            print " ... no match found for %s, appending %s" % ( geoFile, files)
            pass
    
    #for geoFile,prodFile in zip(geoList,prodList):
        #print geoFile,prodFile
    return geoList,prodList

###################################################
#              Granulation Functions              #
###################################################

def _run_crefl(crefl_exe, geoFile, redFile, greenFile, blueFile, outFile):

    currDir = path.abspath('.')
    crefl_exe = path.abspath(path.expanduser(crefl_exe))
    crefl_dir = path.dirname(crefl_exe)

    geoFile = path.abspath(path.expanduser(geoFile))
    blueFile = path.abspath(path.expanduser(blueFile))
    greenFile = path.abspath(path.expanduser(greenFile))
    redFile = path.abspath(path.expanduser(redFile))
    outFile = path.abspath(path.expanduser(outFile))

    LOG.info('outfile = %s' % (outFile))

    os.chdir(crefl_dir)

    try :
        LOG.info('Computing corrected reflectances.')
        cmdStr = '%s viirs --verbose --of=%s %s %s %s %s' % (crefl_exe, outFile, geoFile, blueFile, greenFile, redFile)
        LOG.debug('\t%s' % (cmdStr))
        args = shlex.split(cmdStr)
        LOG.debug('\t%s' % (repr(args)))

        procRetVal = 0
        procObj = subprocess.Popen(args,bufsize=0, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
        procObj.wait()
        procRetVal = procObj.returncode

        procOutput = procObj.stdout.read()

        LOG.debug("%s" % (procOutput))
        
        if not (procRetVal == 0) :
            LOG.error('Computing corrected reflectances failed for crefl command...\n%s.' % (cmdStr))
            sys.exit(procRetVal)

    except Exception, err:
        LOG.warn( "%s" % (str(err)))

    os.chdir(currDir)


def gran_Sdr(geoFile,redFile,greenFile,blueFile,isRadiance=False):
    '''
    Returns the granulated SDR dataset
    '''

    try :
        reload(VIIRS.viirsSDR)
        reload(VIIRS.ViirsData)
        del(viirsSdrObj)
        del(latArr)
        del(lonArr)
        del(redArr)
        del(greenArr)
        del(blueArr)
    except :
        pass

    LOG.info("Creating viirsSdrObj...")
    reload(VIIRS.viirsSDR)
    viirsSdrObj = VIIRS.viirsSDR.viirsSDR()
    LOG.info("done")

    # Determine the correct fillValue
    trimObj = VIIRS.ViirsData.ViirsTrimTable()
    eps = 1.e-6

    rgbKeys = ['lat','lon','lat_0','lon_0','red','green','blue','redName','greenName','blueName']
    RGBfiles = {'red':redFile,'green':greenFile,'blue':blueFile}
    RGBdict = {}

    for chans,chanFileName in zip(['red','green','blue'],['redName','greenName','blueName']):
        LOG.info("Ingesting %s channel ..." % (chans))
        retArr = viirsSdrObj.ingest(geoFile,RGBfiles[chans],isRadiance,1)
        LOG.info("done\n")
        try :

            latArr  = viirsSdrObj.Lat[:,:]
            lonArr  = viirsSdrObj.Lon[:,:]
            ModeGran = viirsSdrObj.ModeGran
            
            lat_0 = latArr[np.shape(latArr)[0]/2,np.shape(latArr)[1]/2]
            lon_0 = lonArr[np.shape(lonArr)[0]/2,np.shape(lonArr)[1]/2]

            RGBdict['lat_0'] = lat_0
            RGBdict['lon_0'] = lon_0

            badGeo = False
            if not (-90. <= lat_0 <= 90.) :
                LOG.error("Latitude of granule midpoint (%f) does not satisfy (-90. <= lat_0 <= 90.)\nfor file %s\n\taborting..." % (lat_0,geoFile))
                badGeo = True
            if not (-180. <= lat_0 <= 180.) :
                LOG.error("Longitude of granule midpoint (%f) does not satisfy (-180. <= lon_0 <= 180.)\nfor file %s\n\taborting..." % (lon_0,geoFile))
                badGeo = True

            if badGeo :
                sys.exit(1)
                
            RGBdict['lat'] = latArr[:,:]
            RGBdict['lon'] = lonArr[:,:]
            RGBdict['ModeGran'] = ModeGran
            RGBdict[chans] = viirsSdrObj.ViirsSDR[:,:]
            RGBdict[chanFileName] = RGBfiles[chans]

        except :
            LOG.error("There was an exception...")
            sys.exit(1)

    return RGBdict


def gran_Sdr2(geoFile,redFile,greenFile,blueFile,isRadiance=False):
    '''
    Returns the granulated SDR dataset
    '''

    try :
        reload(VIIRS.viirsSDR)
        reload(VIIRS.ViirsData)
        del(viirsSdrObj)
        del(latArr)
        del(lonArr)
        del(redArr)
        del(greenArr)
        del(blueArr)
    except :
        pass

    LOG.info("Creating viirsSdrObj...")
    reload(VIIRS.viirsSDR)
    viirsSdrObj = VIIRS.viirsSDR.viirsSDR()
    LOG.info("done")

    # Determine the correct fillValue
    trimObj = VIIRS.ViirsData.ViirsTrimTable()
    eps = 1.e-6

    rgbKeys = ['lat','lon','lat_0','lon_0','red','green','blue','redName','greenName','blueName']
    RGBfiles = {'red':redFile,'green':greenFile,'blue':blueFile}
    RGBdict = {}

    # Ingest the VIIRS geolocation
    LOG.info("Ingesting geolocation...")
    retArr = viirsSdrObj.ingest(geoFile,None,isRadiance,1)

    latArr  = viirsSdrObj.Lat[:,:]
    lonArr  = viirsSdrObj.Lon[:,:]
    ModeGran = viirsSdrObj.ModeGran
    
    lat_0 = latArr[np.shape(latArr)[0]/2,np.shape(latArr)[1]/2]
    lon_0 = lonArr[np.shape(lonArr)[0]/2,np.shape(lonArr)[1]/2]

    RGBdict['lat_0'] = lat_0
    RGBdict['lon_0'] = lon_0

    badGeo = False
    if not (-90. <= lat_0 <= 90.) :
        LOG.error("Latitude of granule midpoint (%f) does not satisfy (-90. <= lat_0 <= 90.)\nfor file %s\n\taborting..." % (lat_0,geoFile))
        badGeo = True
    if not (-180. <= lat_0 <= 180.) :
        LOG.error("Longitude of granule midpoint (%f) does not satisfy (-180. <= lon_0 <= 180.)\nfor file %s\n\taborting..." % (lon_0,geoFile))
        badGeo = True

    if badGeo :
        return 1
        
    RGBdict['lat'] = latArr[:,:]
    RGBdict['lon'] = lonArr[:,:]
    RGBdict['ModeGran'] = ModeGran

    # Ingest the VIIRS red, green and blue channels
    for chans,chanFileName in zip(['red','green','blue'],['redName','greenName','blueName']):
        LOG.info("Ingesting %s channel ..." % (chans))
        retArr = viirsSdrObj.ingest(None,RGBfiles[chans],isRadiance,1)
        LOG.info("done\n")
        try :

            RGBdict[chans] = viirsSdrObj.ViirsSDR[:,:]
            RGBdict[chanFileName] = RGBfiles[chans]

        except :
            LOG.error("There was an exception...")
            return 1

    return RGBdict

def gran_crefl(creflList,geoList):
    '''
    Returns the granulated corrected reflectance
    '''

    try :
        reload(VIIRS.viirsSDR)
        reload(VIIRS.ViirsData)
        del(viirsSdrObj)
        del(latArr)
        del(lonArr)
        del(redArr)
        del(greenArr)
        del(blueArr)
    except :
        pass

    # Determine the correct fillValue
    trimObj = VIIRS.ViirsData.ViirsTrimTable()
    trimMask = trimObj.createModTrimArray(nscans=48,trimType=bool)
    eps = 1.e-6

    rgbKeys = ['lat','lon','lat_0','lon_0','red','green','blue','redName','greenName','blueName']
    creflChanNames = {'red':'CorrRefl_05', 'green':'CorrRefl_04', 'blue':'CorrRefl_03'}

    LOG.info("Creating viirsSdrObj...")
    reload(VIIRS.viirsSDR)
    viirsSdrObj = VIIRS.viirsSDR.viirsSDR()

    #print creflList
    RGBdict = {}
    

    #for grans in np.arange(len(creflList)):
    for grans in np.arange(2):

        creflFile = creflList[grans]
        
        LOG.info("Opening VIIRS crefl file %s..."%(creflFile))
        
        try :
            creflFileObj = SD( creflFile, SDC.READ )
            LOG.info("Successfully opened crefl file")
        except pyHDF.HDF4Error:
            LOG.error("Could not open crefl file: %s" % (creflFile))
            sys.exit(1)

        # Attempt to read geolocation from crefl file. If this fails, 
        # fall back to reading from VIIRS geolocation file

        #LOG.info("Ingesting VIIRS geolocation from crefl file %s..." % (creflFile))
        try :
            latArr  = np.vstack((latArr, creflFileObj.select('Latitude').get()  ))
            lonArr  = np.vstack((lonArr, creflFileObj.select('Longitude').get() ))
            ModeGran = 1
        except NameError :
            latArr = creflFileObj.select('Latitude').get()
            lonArr = creflFileObj.select('Longitude').get()
            ModeGran = 1
        
        LOG.info("Intermediate latArr.shape = %s" % (str(latArr.shape)))
        LOG.info("Intermediate lonArr.shape = %s" % (str(lonArr.shape)))

        for chans in ['red','green','blue']:
            colourArr = '%sArr'%(chans)
            LOG.info("Processing %s..." % (colourArr))
            try :
                chanArr = creflFileObj.select(creflChanNames[chans]).get()
                LOG.info("%s chanArr shape is %s" % (chans,str(chanArr.shape)))
                vars()[colourArr] = np.vstack((vars()[colourArr],chanArr))
                LOG.info("Finished reading VIIRS corrected reflectance %s from %s" % (chans,creflFile))
            except KeyError:
                LOG.info("Reading first VIIRS corrected reflectance %s from %s..." % (creflChanNames[chans],creflFile))
                vars()[colourArr]  = creflFileObj.select(creflChanNames[chans]).get()
                LOG.info("Finished reading first VIIRS corrected reflectance %s from %s" % (chans,creflFile))

            #LOG.info("Intermediate %sArr.shape = %s" % (chans,str(vars()[chans+'Arr'].shape)))

        #try :
        #except Exception, err:
            #creflFileObj.end()
            #LOG.error( "%s" % (str(err)))
            #LOG.error( "Unable to read VIIRS corrected reflectance %s from %s, aborting..." % (chans,creflFile))
            #sys.exit(1)


        # Close crefl file
        creflFileObj.end()
        print ""

    print redArr

    #except Exception, err:
        #LOG.error( "%s" % (str(err)))
    
    lat_0 = latArr[np.shape(latArr)[0]/2,np.shape(latArr)[1]/2]
    lon_0 = lonArr[np.shape(lonArr)[0]/2,np.shape(lonArr)[1]/2]

    RGBdict['lat_0'] = lat_0
    RGBdict['lon_0'] = lon_0
    RGBdict['lat'] = latArr[:,:]
    RGBdict['lon'] = lonArr[:,:]
    RGBdict['ModeGran'] = ModeGran
    RGBdict['redName'] = 'SVM05_npp_d20120705_t2028301_e2029543_b03568_c20120706035000918229_noaa_ops.h5'
    RGBdict['greenName'] = 'SVM04_npp_d20120705_t2028301_e2029543_b03568_c20120706034958022764_noaa_ops.h5'
    RGBdict['blueName'] = 'SVM03_npp_d20120705_t2028301_e2029543_b03568_c20120706035000894469_noaa_ops.h5'

    # TODO : Impose proper bow-tie deletion on crefl output.


    return RGBdict


def _create_map_image_RGB(RGBdict,pointSize=1.,scale=1.3,mapRes='c',stride=1,\
    sdrKey='M15',isRadiance=False,prodFileName='',outFileName='out.png',dpi=300,titleStr='VIIRS RGB'):
    '''
    Plots the VIIRS Scientific Data Record (SDR)on an orthographic projection
    '''

    gridLat = RGBdict['lat']
    gridLon = RGBdict['lon']
    lat_0 = RGBdict['lat_0']
    lon_0 = RGBdict['lon_0']
    ModeGran = RGBdict['ModeGran']

    # Setup plotting data
    reload(VIIRS.ViirsData)

    # If we have a zero size array, make a dummy dataset
    # to span the allowed data range, which will be plotted with 
    # vanishing pointsize
    if (np.shape(gridLon)[0]==0) :
        print "We have no valid data, synthesising dummy data..."
        gridLat = np.array([lat_0,lat_0])
        gridLon = np.array([lon_0,lon_0])
        gridData = np.array([vmin,vmax])
        pointSize = 0.001

    figWidth = 5. # inches
    figHeight = 4. # inches

    # Create figure with default size, and create canvas to draw on
    fig = Figure(figsize=((figWidth,figHeight)))
    canvas = FigureCanvas(fig)

    # Create main axes instance, leaving room for colorbar at bottom,
    # and also get the Bbox of the axes instance
    ax_rect = [0.05, 0.18, 0.9, 0.75  ] # [left,bottom,width,height]
    ax = fig.add_axes(ax_rect)

    # Granule axis title
    ax_title = ppl.setp(ax,title=prodFileName)
    ppl.setp(ax_title,fontsize=6)
    ppl.setp(ax_title,family="monospace")

    # Create Basemap instance
    # set 'ax' keyword so pylab won't be imported.
    #m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,ax=ax,fix_aspect=True,resolution='c')
    print "scale = ",scale
    m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,\
        ax=ax,fix_aspect=True,resolution=mapRes,\
        llcrnrx = -1. * scale * 3200. * 750./2.,\
        llcrnry = -1. * scale * 3200. * 750./2.,\
        urcrnrx =       scale * 3200. * 750./2.,\
        urcrnry =       scale * 3200. * 750./2.)

    x,y=m(np.array(gridLon),np.array(gridLat))

    # Some map style configuration stufff
    #m.drawlsmask(ax=ax,land_color='gray',ocean_color='black',lakes=True)
    m.drawmapboundary(ax=ax,linewidth=0.01,fill_color='grey')
    m.drawcoastlines(ax=ax,linewidth=0.3,color='white')
    #m.fillcontinents(ax=ax,color='gray',lake_color='black',zorder=0)
    #m.drawparallels(np.arange(-90.,120.,30.),color='white')
    #m.drawmeridians(np.arange(0.,420.,60.),color='white')

    if (ModeGran == 0) :
        fig.text(0.5, 0.555, 'NIGHT',fontsize=30, color='white',ha='center',va='center',alpha=0.6)
    else :

        # Convert the scaledImg array to something that can be displayed using scatter...

        scaledImg = RGBdict['scaledImg']

        trimObj = VIIRS.ViirsData.ViirsTrimTable()
        trimMask = trimObj.createModTrimArray(nscans=48,trimType=bool)
        raveldTrimMask = np.ravel(trimMask)
        onboard_pt_value = trimObj.sdrTypeFill['ONBOARD_PT_FILL'][scaledImg.dtype.name]

        # Mask the geolocation
        x = ma.array(np.ravel(x),mask=raveldTrimMask)
        y = ma.array(np.ravel(y),mask=raveldTrimMask)
        
        # Scale channels to be 0..1
        R = scaledImg[:,:,0].astype('float')/255.
        G = scaledImg[:,:,1].astype('float')/255.
        B = scaledImg[:,:,2].astype('float')/255.

        # Ravel the channel arrays
        R = np.ravel(R)
        G = np.ravel(G)
        B = np.ravel(B)

        # Mask the RGB arrays
        R = ma.array(R,mask=raveldTrimMask)
        G = ma.array(G,mask=raveldTrimMask)
        B = ma.array(B,mask=raveldTrimMask)

        # Compress the geolocation and RGB arrays
        x = x.compressed()
        y = y.compressed()
        R = R.compressed()
        G = G.compressed()
        B = B.compressed()

        # Shrink the geolocation and RGB arrays
        print "Stride = ",stride
        print "dpi = ",dpi
        x = x[::stride]
        y = y[::stride]
        R = R[::stride]
        G = G[::stride]
        B = B[::stride]

        # Construct the data array
        gridData = np.transpose(np.array([R,G,B]))

        print x.shape
        print y.shape
        print gridData.shape

        # Plot the granule data
        # TODO: Use pcolormesh instead.
        cs = m.scatter(x,y,s=pointSize,c=gridData,axes=ax,edgecolors='none',vmin=0.,vmax=1.)

    # RGB information at the bottom of the plot.

    timeString = 'Creation date: %s' %(datetime.strftime(datetime.utcnow(),"%Y-%m-%d %H:%M:%S Z"))
    fig.text(0.98, 0.01, timeString,fontsize=5, color='gray',ha='right',va='bottom',alpha=0.9)

    channelDict = {'red':string.split(path.basename(RGBdict['redName']),'_')[0][2:],\
                   'green':string.split(path.basename(RGBdict['greenName']),'_')[0][2:],\
                   'blue':string.split(path.basename(RGBdict['blueName']),'_')[0][2:]}

    RGBstring = "R : %s (%s $\mu\mathrm{m}$); G : %s (%s $\mu\mathrm{m}$); B : %s (%s $\mu\mathrm{m}$)" % \
            (channelDict['red'],\
             str(VIIRS.ViirsData.ViirsBandCenter[channelDict['red']]),\
             channelDict['green'],\
             str(VIIRS.ViirsData.ViirsBandCenter[channelDict['green']]),\
             channelDict['blue'],\
             str(VIIRS.ViirsData.ViirsBandCenter[channelDict['blue']])
            )
    fig.text(0.5, 0.07, RGBstring,fontsize=10, color='black',ha='center',va='bottom',alpha=1.0)


    #
    # Add a small globe with the swath indicated on it #
    #

    # Create main axes instance, leaving room for colorbar at bottom,
    # and also get the Bbox of the axes instance
    glax_rect = [0.81, 0.75, 0.18, 0.20 ] # [left,bottom,width,height]
    glax = fig.add_axes(glax_rect)

    m_globe = Basemap(lat_0=0.,lon_0=0.,\
        ax=glax,resolution='c',area_thresh=10000.,projection='robin')

    # If we previously had a zero size data array, increase the pointSize
    # so the data points are visible on the global plot
    pointSize = 0.2
    if (np.shape(gridLon)[0]==2) :
        print "We have no valid data, synthesising dummy data..."
        pointSize = 5.

    x,y = m_globe(np.array(gridLon),np.array(gridLat))
    swath = np.zeros(np.shape(x),dtype=int)

    m_globe.drawcoastlines(ax=glax,linewidth=0.1)
    m_globe.fillcontinents(ax=glax,color='gray',zorder=0)
    m_globe.drawmapboundary(linewidth=0.1)

    p_globe = m_globe.scatter(x,y,s=pointSize,c="red",axes=glax,edgecolors='none')

    # Globe axis title
    glax_xlabel = ppl.setp(glax,xlabel=titleStr)
    ppl.setp(glax_xlabel,fontsize=6)

    # Redraw the figure
    canvas.draw()

    # save image 
    print "Writing file to ",outFileName
    canvas.print_figure(outFileName,dpi=dpi)


def _create_map_image_RGB_noravel(RGBdict,pointSize=1.,scale=1.3,mapRes='c',stride=1,\
    sdrKey='M15',isRadiance=False,prodFileName='',outFileName='out.png',dpi=300,titleStr='VIIRS RGB'):
    '''
    Plots the VIIRS Scientific Data Record (SDR)on an orthographic projection
    '''

    gridLat = RGBdict['lat']
    gridLon = RGBdict['lon']
    lat_0 = RGBdict['lat_0']
    lon_0 = RGBdict['lon_0']
    ModeGran = RGBdict['ModeGran']

    # Setup plotting data
    reload(VIIRS.ViirsData)

    # If we have a zero size array, make a dummy dataset
    # to span the allowed data range, which will be plotted with 
    # vanishing pointsize
    if (np.shape(gridLon)[0]==0) :
        print "We have no valid data, synthesising dummy data..."
        gridLat = np.array([lat_0,lat_0])
        gridLon = np.array([lon_0,lon_0])
        gridData = np.array([vmin,vmax])
        pointSize = 0.001

    figWidth = 5. # inches
    figHeight = 4. # inches

    # Create figure with default size, and create canvas to draw on
    fig = Figure(figsize=((figWidth,figHeight)))
    canvas = FigureCanvas(fig)

    # Create main axes instance, leaving room for colorbar at bottom,
    # and also get the Bbox of the axes instance
    ax_rect = [0.05, 0.18, 0.9, 0.75  ] # [left,bottom,width,height]
    ax = fig.add_axes(ax_rect)

    # Granule axis title
    ax_title = ppl.setp(ax,title=prodFileName)
    ppl.setp(ax_title,fontsize=6)
    ppl.setp(ax_title,family="monospace")

    # Create Basemap instance
    # set 'ax' keyword so pylab won't be imported.
    #m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,ax=ax,fix_aspect=True,resolution='c')
    print "scale = ",scale
    m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,\
        ax=ax,fix_aspect=True,resolution=mapRes,\
        llcrnrx = -1. * scale * 3200. * 750./2.,\
        llcrnry = -1. * scale * 3200. * 750./2.,\
        urcrnrx =       scale * 3200. * 750./2.,\
        urcrnry =       scale * 3200. * 750./2.)

    x,y=m(np.array(gridLon),np.array(gridLat))

    # Some map style configuration stufff
    #m.drawlsmask(ax=ax,land_color='gray',ocean_color='black',lakes=True)
    m.drawmapboundary(ax=ax,linewidth=0.01,fill_color='grey')
    m.drawcoastlines(ax=ax,linewidth=0.3,color='white')
    #m.fillcontinents(ax=ax,color='gray',lake_color='black',zorder=0)
    #m.drawparallels(np.arange(-90.,120.,30.),color='white')
    #m.drawmeridians(np.arange(0.,420.,60.),color='white')

    if (ModeGran == 0) :
        fig.text(0.5, 0.555, 'NIGHT',fontsize=30, color='white',ha='center',va='center',alpha=0.6)
    else :

        # Convert the scaledImg array to something that can be displayed using scatter...

        scaledImg = RGBdict['scaledImg']

        trimObj = VIIRS.ViirsData.ViirsTrimTable()
        trimMask = trimObj.createModTrimArray(nscans=48,trimType=bool)
        raveldTrimMask = np.ravel(trimMask)
        onboard_pt_value = trimObj.sdrTypeFill['ONBOARD_PT_FILL'][scaledImg.dtype.name]

        # Mask the geolocation
        x = ma.array(np.ravel(x),mask=raveldTrimMask)
        y = ma.array(np.ravel(y),mask=raveldTrimMask)
        
        # Scale channels to be 0..1
        R = scaledImg[:,:,0].astype('float')/255.
        G = scaledImg[:,:,1].astype('float')/255.
        B = scaledImg[:,:,2].astype('float')/255.

        # Ravel the channel arrays
        R = np.ravel(R)
        G = np.ravel(G)
        B = np.ravel(B)

        # Mask the RGB arrays
        R = ma.array(R,mask=raveldTrimMask)
        G = ma.array(G,mask=raveldTrimMask)
        B = ma.array(B,mask=raveldTrimMask)

        # Compress the geolocation and RGB arrays
        x = x.compressed()
        y = y.compressed()
        R = R.compressed()
        G = G.compressed()
        B = B.compressed()

        # Shrink the geolocation and RGB arrays
        print "Stride = ",stride
        print "dpi = ",dpi
        x = x[::stride]
        y = y[::stride]
        R = R[::stride]
        G = G[::stride]
        B = B[::stride]

        # Construct the data array
        gridData = np.transpose(np.array([R,G,B]))

        print x.shape
        print y.shape
        print gridData.shape

        # Plot the granule data
        # TODO: Use pcolormesh instead.
        cs = m.scatter(x,y,s=pointSize,c=gridData,axes=ax,edgecolors='none',vmin=0.,vmax=1.)

    # RGB information at the bottom of the plot.

    timeString = 'Creation date: %s' %(datetime.strftime(datetime.utcnow(),"%Y-%m-%d %H:%M:%S Z"))
    fig.text(0.98, 0.01, timeString,fontsize=5, color='gray',ha='right',va='bottom',alpha=0.9)

    channelDict = {'red':string.split(path.basename(RGBdict['redName']),'_')[0][2:],\
                   'green':string.split(path.basename(RGBdict['greenName']),'_')[0][2:],\
                   'blue':string.split(path.basename(RGBdict['blueName']),'_')[0][2:]}

    RGBstring = "R : %s (%s $\mu\mathrm{m}$); G : %s (%s $\mu\mathrm{m}$); B : %s (%s $\mu\mathrm{m}$)" % \
            (channelDict['red'],\
             str(VIIRS.ViirsData.ViirsBandCenter[channelDict['red']]),\
             channelDict['green'],\
             str(VIIRS.ViirsData.ViirsBandCenter[channelDict['green']]),\
             channelDict['blue'],\
             str(VIIRS.ViirsData.ViirsBandCenter[channelDict['blue']])
            )
    fig.text(0.5, 0.07, RGBstring,fontsize=10, color='black',ha='center',va='bottom',alpha=1.0)


    #
    # Add a small globe with the swath indicated on it #
    #

    # Create main axes instance, leaving room for colorbar at bottom,
    # and also get the Bbox of the axes instance
    glax_rect = [0.81, 0.75, 0.18, 0.20 ] # [left,bottom,width,height]
    glax = fig.add_axes(glax_rect)

    m_globe = Basemap(lat_0=0.,lon_0=0.,\
        ax=glax,resolution='c',area_thresh=10000.,projection='robin')

    # If we previously had a zero size data array, increase the pointSize
    # so the data points are visible on the global plot
    pointSize = 0.2
    if (np.shape(gridLon)[0]==2) :
        print "We have no valid data, synthesising dummy data..."
        pointSize = 5.

    x,y = m_globe(np.array(gridLon),np.array(gridLat))
    swath = np.zeros(np.shape(x),dtype=int)

    m_globe.drawcoastlines(ax=glax,linewidth=0.1)
    m_globe.fillcontinents(ax=glax,color='gray',zorder=0)
    m_globe.drawmapboundary(linewidth=0.1)

    p_globe = m_globe.scatter(x,y,s=pointSize,c="red",axes=glax,edgecolors='none')

    # Globe axis title
    glax_xlabel = ppl.setp(glax,xlabel=titleStr)
    ppl.setp(glax_xlabel,fontsize=6)

    # Redraw the figure
    canvas.draw()

    # save image 
    print "Writing file to ",outFileName
    canvas.print_figure(outFileName,dpi=dpi)


class BandData():
    def __init__(self,geoGlob):

        geoFile = path.abspath(glob(geoGlob)[0])
        dateStamp = string.split(path.basename(geoFile),'_')[2]
        timeStamp = string.split(path.basename(geoFile),'_')[3]

        print geoFile
        print dateStamp
        print timeStamp

        bands = ['03','04','05','10','11']

        bandFiles = {}
        for band in bands:
            bandGlob = string.replace(geoGlob,'GMODO','SVM'+band)
            bandFiles['SVM'+band] = path.abspath(glob(bandGlob)[0])
            print bandFiles['SVM'+band]


def viirs_rgb(imgType,RGBdict,swathImg,swathName=None):

    print 'lat_0 = %f' % (RGBdict['lat_0'])
    print 'lon_0 = %f' % (RGBdict['lon_0'])
    print 'lat.shape = %s' % (repr(RGBdict['lat'].shape))
    print 'lon.shape = %s' % (repr(RGBdict['lon'].shape))
    print 'red.shape = %s' % (repr(RGBdict['red'].shape))
    print 'green.shape = %s' % (repr(RGBdict['green'].shape))
    print 'blue.shape = %s' % (repr(RGBdict['blue'].shape))

    # Construct the Bow-tie deletion mask
    trimObj = VIIRS.ViirsData.ViirsTrimTable()
    trimMask = trimObj.createModTrimArray(nscans=48,trimType=bool)

    R = ma.masked_less(RGBdict['red'],-1.5)
    G = ma.masked_less(RGBdict['green'],-1.5)
    B = ma.masked_less(RGBdict['blue'],-1.5)

    print "min RGBdict['red'] = %f" % (np.min(R))
    print "min RGBdict['green'] = %f" % (np.min(G))
    print "min RGBdict['blue'] = %f" % (np.min(B))
    print "max RGBdict['red'] = %f" % (np.max(R))
    print "max RGBdict['green'] = %f" % (np.max(G))
    print "max RGBdict['blue'] = %f" % (np.max(B))

    bandMin = max(np.min(R),np.min(G),np.min(B))
    bandMax = max(np.max(R),np.max(G),np.max(B))

    band_R_scaled = 255. * (R - np.min(R)) / np.max(R - np.min(R))
    band_G_scaled = 255. * (G - np.min(G)) / np.max(G - np.min(G))
    band_B_scaled = 255. * (B - np.min(B)) / np.max(B - np.min(B))

    newimg = np.zeros((768,3200,3),dtype=np.float32)
    newimg[:,:,0] = band_R_scaled[:,:]
    newimg[:,:,1] = band_G_scaled[:,:]
    newimg[:,:,2] = band_B_scaled[:,:]

    # Make an image object from a numpy array
    im = Image.fromarray(newimg.astype(np.uint8))

    # Limits for differential scaling

    if imgType == 'TrueColor' or imgType == 'Vegetation' :
        # True Colour
        input  = np.array([0,  30,  60, 120, 190, 255],dtype=np.float32)
        #output = np.array([0, 110, 160, 210, 240, 255],dtype=np.float32) # Original scaling
        output = np.array([0, 180, 200, 220, 240, 255],dtype=np.float32) # New scaling
    elif imgType == 'SnowIce' :
        # Snow/Ice
        input  = np.array([0, 25,  55, 100, 255],dtype=np.float32)
        output = np.array([0, 90, 140, 175, 255],dtype=np.float32)
    elif imgType == 'CloudPhase' :
        newimg[:,:,1] = 255. - band_B_scaled
        input  = np.array([0,  30,  60, 120, 190, 255],dtype=np.float32)
        output = np.array([0,  30,  60, 120, 190, 255],dtype=np.float32) # New scaling
    else :
        print "error: invalid band combination. Try 'TrueColor','SnowIce','Vegetation'"
        return 1

    inputRange = input[1 :]-input[:-1]
    outputRange = output[1 :]-output[:-1]
    inputMin,inputMax = input[:-1], input[1 :]
    outputMin,outputMax = output[:-1], output[1 :]

    slope = (outputMax-outputMin)/(inputMax-inputMin)
    intercept = outputMax - (slope*inputMax)

    scaledImg = np.zeros((768,3200,3),dtype=np.float32)

    print "Scaling red channel..."
    for pair in np.arange(np.size(inputRange)):
        inputIdx = np.where((newimg[:,:,0]>=inputMin[pair])*(newimg[:,:,0]<inputMax[pair]))
        rows,cols = inputIdx[0],inputIdx[1]
        print "min,max of unscaled red channel:",np.min(newimg[rows,cols,0]),np.max(newimg[rows,cols,0])
        scaledImg[rows,cols,0] = slope[pair] * newimg[rows,cols,0] + intercept[pair]
        print "min,max of scaled red channel:",np.min(newimg[rows,cols,0]),np.max(newimg[rows,cols,0])

    print "\nScaling green channel..."
    for pair in np.arange(np.size(inputRange)):
        inputIdx = np.where((newimg[:,:,1]>=inputMin[pair])*(newimg[:,:,1]<inputMax[pair]))
        rows,cols = inputIdx[0],inputIdx[1]
        print "min,max of unscaled green channel:",np.min(newimg[rows,cols,1]),np.max(newimg[rows,cols,1])
        scaledImg[rows,cols,1] = slope[pair] * newimg[rows,cols,1] + intercept[pair]
        print "min,max of scaled green channel:",np.min(newimg[rows,cols,1]),np.max(newimg[rows,cols,1])

    print "\nScaling blue channel..."
    for pair in np.arange(np.size(inputRange)):
        inputIdx = np.where((newimg[:,:,2]>=inputMin[pair])*(newimg[:,:,2]<inputMax[pair]))
        rows,cols = inputIdx[0],inputIdx[1]
        print "min,max of unscaled blue channel:",np.min(newimg[rows,cols,2]),np.max(newimg[rows,cols,2])
        scaledImg[rows,cols,2] = slope[pair] * newimg[rows,cols,2] + intercept[pair]
        print "min,max of scaled blue channel:",np.min(newimg[rows,cols,2]),np.max(newimg[rows,cols,2])

    scaledImg = scaledImg.astype('uint8')

    RGB_trimMask = np.zeros((768,3200,3),dtype=np.bool)
    RGB_trimMask[:,:,0] = trimMask[:,:]
    RGB_trimMask[:,:,1] = trimMask[:,:]
    RGB_trimMask[:,:,2] = trimMask[:,:]

    onboard_pt_value = trimObj.sdrTypeFill['ONBOARD_PT_FILL']['uint8']
    scaledImg_masked = np.zeros((768,3200,3),dtype=np.uint8)
    scaledImg_masked = ma.array(scaledImg,mask=RGB_trimMask,fill_value=onboard_pt_value)

    # Make some pixels red...
    afires_mask = np.zeros((768,3200,3),dtype=np.uint8)
    afires_mask[:,:,0] = 255 * np.ones((768,3200),dtype=np.uint8)

    scaledImg_masked[380:390,1600:1610,:] = afires_mask[380:390,1600:1610,:] 
    scaledImg_masked[8,2264,:] = afires_mask[8,2264,:] 
    scaledImg_masked[9,2491,:] = afires_mask[9,2491,:] 

    if swathImg :
        im = Image.fromarray(scaledImg_masked[::-1,::-1,:])
        LOG.info("Initial swathName = %r" % (swathName))
        if swathName == None :
            swathName = 'Viirs_swath_truecolor.png'
        else :
            swathName = string.replace(swathName,'.png','_swath.png')
        LOG.info("Final swathName = %r" % (swathName))
        im.save(swathName)

    scaledImg = scaledImg_masked.filled()

    RGBdict['scaledImg'] = scaledImg

    return RGBdict


def crefl_rgb(imgType,RGBdict,swathImg,swathName=None):

    print 'lat_0 = %f' % (RGBdict['lat_0'])
    print 'lon_0 = %f' % (RGBdict['lon_0'])
    print 'lat.shape = %s' % (repr(RGBdict['lat'].shape))
    print 'lon.shape = %s' % (repr(RGBdict['lon'].shape))
    print 'red.shape = %s' % (repr(RGBdict['red'].shape))
    print 'green.shape = %s' % (repr(RGBdict['green'].shape))
    print 'blue.shape = %s' % (repr(RGBdict['blue'].shape))

    # Construct the Bow-tie deletion mask
    trimObj = VIIRS.ViirsData.ViirsTrimTable()
    trimMask = trimObj.createModTrimArray(nscans=48,trimType=bool)

    R = ma.masked_equal(RGBdict['red'],32767)
    G = ma.masked_equal(RGBdict['green'],32767)
    B = ma.masked_equal(RGBdict['blue'],32767)

    print "min RGBdict['red'] = %f" % (np.min(R))
    print "min RGBdict['green'] = %f" % (np.min(G))
    print "min RGBdict['blue'] = %f" % (np.min(B))
    print "max RGBdict['red'] = %f" % (np.max(R))
    print "max RGBdict['green'] = %f" % (np.max(G))
    print "max RGBdict['blue'] = %f" % (np.max(B))

    bandMin = max(np.min(R),np.min(G),np.min(B))
    bandMax = max(np.max(R),np.max(G),np.max(B))

    band_R_scaled = 255. * (R - np.min(R)) / np.max(R - np.min(R))
    band_G_scaled = 255. * (G - np.min(G)) / np.max(G - np.min(G))
    band_B_scaled = 255. * (B - np.min(B)) / np.max(B - np.min(B))

    newimg = np.zeros((768,3200,3),dtype=np.float32)
    newimg[:,:,0] = band_R_scaled[:,:]
    newimg[:,:,1] = band_G_scaled[:,:]
    newimg[:,:,2] = band_B_scaled[:,:]

    # Make an image object from a numpy array
    im = Image.fromarray(newimg.astype(np.uint8))

    # Limits for differential scaling

    if imgType == 'TrueColor' or imgType == 'Vegetation' :
        # True Colour
        input  = np.array([0,  30,  60, 120, 190, 255],dtype=np.float32)
        #output = np.array([0, 110, 160, 210, 240, 255],dtype=np.float32) # Original scaling
        output = np.array([0, 180, 200, 220, 240, 255],dtype=np.float32) # New scaling
    elif imgType == 'SnowIce' :
        # Snow/Ice
        input  = np.array([0, 25,  55, 100, 255],dtype=np.float32)
        output = np.array([0, 90, 140, 175, 255],dtype=np.float32)
    elif imgType == 'CloudPhase' :
        newimg[:,:,1] = 255. - band_B_scaled
        input  = np.array([0,  30,  60, 120, 190, 255],dtype=np.float32)
        output = np.array([0,  30,  60, 120, 190, 255],dtype=np.float32) # New scaling
    else :
        print "error: invalid band combination. Try 'TrueColor','SnowIce','Vegetation'"
        return 1

    inputRange = input[1 :]-input[:-1]
    outputRange = output[1 :]-output[:-1]
    inputMin,inputMax = input[:-1], input[1 :]
    outputMin,outputMax = output[:-1], output[1 :]

    slope = (outputMax-outputMin)/(inputMax-inputMin)
    intercept = outputMax - (slope*inputMax)

    scaledImg = np.zeros((768,3200,3),dtype=np.float32)

    print "Scaling red channel..."
    for pair in np.arange(np.size(inputRange)):
        inputIdx = np.where((newimg[:,:,0]>=inputMin[pair])*(newimg[:,:,0]<inputMax[pair]))
        rows,cols = inputIdx[0],inputIdx[1]
        print "min,max of unscaled red channel:",np.min(newimg[rows,cols,0]),np.max(newimg[rows,cols,0])
        scaledImg[rows,cols,0] = slope[pair] * newimg[rows,cols,0] + intercept[pair]
        print "min,max of scaled red channel:",np.min(newimg[rows,cols,0]),np.max(newimg[rows,cols,0])

    print "\nScaling green channel..."
    for pair in np.arange(np.size(inputRange)):
        inputIdx = np.where((newimg[:,:,1]>=inputMin[pair])*(newimg[:,:,1]<inputMax[pair]))
        rows,cols = inputIdx[0],inputIdx[1]
        print "min,max of unscaled green channel:",np.min(newimg[rows,cols,1]),np.max(newimg[rows,cols,1])
        scaledImg[rows,cols,1] = slope[pair] * newimg[rows,cols,1] + intercept[pair]
        print "min,max of scaled green channel:",np.min(newimg[rows,cols,1]),np.max(newimg[rows,cols,1])

    print "\nScaling blue channel..."
    for pair in np.arange(np.size(inputRange)):
        inputIdx = np.where((newimg[:,:,2]>=inputMin[pair])*(newimg[:,:,2]<inputMax[pair]))
        rows,cols = inputIdx[0],inputIdx[1]
        print "min,max of unscaled blue channel:",np.min(newimg[rows,cols,2]),np.max(newimg[rows,cols,2])
        scaledImg[rows,cols,2] = slope[pair] * newimg[rows,cols,2] + intercept[pair]
        print "min,max of scaled blue channel:",np.min(newimg[rows,cols,2]),np.max(newimg[rows,cols,2])

    scaledImg = scaledImg.astype('uint8')

    # TODO : Actually apply the mask to the arrays.

    RGB_trimMask = np.zeros((768,3200,3),dtype=np.bool)
    RGB_trimMask[:,:,0] = trimMask[:,:]
    RGB_trimMask[:,:,1] = trimMask[:,:]
    RGB_trimMask[:,:,2] = trimMask[:,:]

    onboard_pt_value = trimObj.sdrTypeFill['ONBOARD_PT_FILL']['uint8']
    scaledImg_masked = np.zeros((768,3200,3),dtype=np.uint8)
    scaledImg_masked = ma.array(scaledImg,mask=RGB_trimMask,fill_value=onboard_pt_value)

    if swathImg :
        im = Image.fromarray(scaledImg_masked)
        LOG.info("Initial swathName = %r" % (swathName))
        if swathName == None :
            swathName = 'Viirs_swath_truecolor_crefl.png'
        else :
            swathName = string.replace(swathName,'.png','_swath.png')
        LOG.info("Final swathName = %r" % (swathName))
        im.save(swathName)

    scaledImg = scaledImg_masked.filled()

    RGBdict['scaledImg'] = scaledImg

    return RGBdict

######################################################
######################################################
######################################################

def basemap_rgb(RGBdict):

    # it looks like the 
    reload(VIIRS.ViirsData)

    figWidth = 5. # inches
    figHeight = 4. # inches

    # Create figure with default size, and create canvas to draw on
    fig = Figure(figsize=((figWidth,figHeight)))
    canvas = FigureCanvas(fig)

    # Create main axes instance, leaving room for colorbar at bottom,
    # and also get the Bbox of the axes instance
    ax_rect = [0.05, 0.18, 0.9, 0.75  ] # [left,bottom,width,height]
    ax = fig.add_axes(ax_rect)

    # Granule axis title
    ax_title = ppl.setp(ax,title="VIIRS RGB")
    ppl.setp(ax_title,fontsize=6)
    ppl.setp(ax_title,family="monospace")

    # Create Basemap instance
    # set 'ax' keyword so pylab won't be imported.
    #m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,ax=ax,fix_aspect=True,resolution='c')
    scale = 1.3
    lat_0 = RGBdict['lat_0']
    lon_0 = RGBdict['lon_0']

    m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,\
        ax=ax,fix_aspect=True,resolution='l',\
        llcrnrx = -1. * scale * 3200. * 750./2.,\
        llcrnry = -1. * scale * 3200. * 750./2.,\
        urcrnrx =       scale * 3200. * 750./2.,\
        urcrnry =       scale * 3200. * 750./2.)

    x,y=m(np.array(RGBdict['lon']),np.array(RGBdict['lat']))

    # Some map style configuration stufff
    m.drawmapboundary(ax=ax,linewidth=0.01,fill_color='grey')
    m.drawcoastlines(ax=ax,linewidth=0.3,color='white')
    #m.fillcontinents(ax=ax,color='gray',lake_color='black',zorder=0)
    #m.drawparallels(np.arange(-90.,120.,30.),color='white')
    #m.drawmeridians(np.arange(0.,420.,60.),color='white')

    # Convert the scaledImg array to something that can be displayed using scatter...

    scaledImg = RGBdict['scaledImg']

    #shrink = 1
    shrink = 13

    trimObj = VIIRS.ViirsData.ViirsTrimTable()
    trimMask = trimObj.createModTrimArray(nscans=48,trimType=bool)
    raveldTrimMask = np.ravel(trimMask)
    onboard_pt_value = trimObj.sdrTypeFill['ONBOARD_PT_FILL'][scaledImg.dtype.name]

    # Mask the geolocation
    x = ma.array(np.ravel(x),mask=raveldTrimMask)
    y = ma.array(np.ravel(y),mask=raveldTrimMask)
    
    # Scale channels to be 0..1
    R = scaledImg[:,:,0].astype('float')/255.
    G = scaledImg[:,:,1].astype('float')/255.
    B = scaledImg[:,:,2].astype('float')/255.

    # Ravel the channel arrays
    R = np.ravel(R)
    G = np.ravel(G)
    B = np.ravel(B)

    # Mask the RGB arrays
    R = ma.array(R,mask=raveldTrimMask)
    G = ma.array(G,mask=raveldTrimMask)
    B = ma.array(B,mask=raveldTrimMask)

    # Compress the geolocation and RGB arrays
    x = x.compressed()
    y = y.compressed()
    R = R.compressed()
    G = G.compressed()
    B = B.compressed()

    # Shrink the geolocation and RGB arrays
    x = x[::shrink]
    y = y[::shrink]
    R = R[::shrink]
    G = G[::shrink]
    B = B[::shrink]

    # Construct the data array
    data = np.transpose(np.array([R,G,B]))

    print x.shape
    print y.shape
    print data.shape

    pointSize=0.2

    # Plot the granule data
    cs = m.scatter(x,y,s=pointSize,c=data,axes=ax,faceted=False,vmin=0.,vmax=1.)

    # Redraw the figure
    canvas.draw()

    # save image 
    canvas.print_figure("test.png",dpi=300)

###################################################
#                  Main Function                  #
###################################################

def main():

    prodChoices=['TrueColor','SnowIce','CloudPhase','Vegetation']
    mapRes = ['c','l','i']

    description = \
    '''
    This is a brief description of %prog
    '''
    usage = "usage: %prog [mandatory args] [options]"
    version = version="%prog PEATE"
    parser = optparse.OptionParser(description=description,usage=usage,version=version)

    # Mandatory arguments
    mandatoryGroup = optparse.OptionGroup(parser, "Mandatory Arguments",
                        "At a minimum these arguments must be specified")

    parser.add_option_group(mandatoryGroup)

    # Optional arguments
    optionalGroup = optparse.OptionGroup(parser, "Extra Options",
                        "These options may be used to customize plot characteristics.")

    optionalGroup.add_option('-g','--geo_file',
                      action="store",
                      dest="geoFile" ,
                      type="string",
                      help="The full path of the VIIRS geolocation file")
    optionalGroup.add_option('--red_file',
                      action="store",
                      dest="redFile",
                      type="string",
                      help="The full path of the VIIRS SDR file corresponding to the red channel.")
    optionalGroup.add_option('--green_file',
                      action="store",
                      dest="greenFile",
                      type="string",
                      help="The full path of the VIIRS SDR file corresponding to the green channel.")
    optionalGroup.add_option('--blue_file',
                      action="store",
                      dest="blueFile",
                      type="string",
                      help="The full path of the VIIRS SDR file corresponding to the blue channel.")
    optionalGroup.add_option('-r','--svn_revision',
                      action="store",
                      dest="svnRevision",
                      default=string.split(__version__," ")[2],
                      type="string",
                      help="The Subversion revision number/tag of this script")
    optionalGroup.add_option('-R','--svn_repo_path',
                      action="store",
                      dest="svnRepoPath",
                      default="https://svn.ssec.wisc.edu/repos/geoffc/Python/VIIRS/"+path.basename(sys.argv[0]),
                      type="string",
                      help="The full Subversion repository path of this script [default is %default].")
    optionalGroup.add_option('-d','--dpi',
                      action="store",
                      dest="dpi",
                      default='200.',
                      type="float",
                      help="The resolution in dots per inch of the output png file. [default: %default]")
    optionalGroup.add_option('-s','--scale',
                      action="store",
                      dest="scale",
                      default='1.3',
                      type="float",
                      help="The extent of the plot viewport as a proportion of the VIIRS swath width (2400 km). [default: %default]")
    optionalGroup.add_option('-S','--stride',
                      action="store",
                      dest="stride",
                      default='1',
                      type="int",
                      help="Sample every STRIDE pixels in the VIIRS IP/SDR product. [default: %default]")
    optionalGroup.add_option('-P','--point_size',
                      action="store",
                      dest="pointSize",
                      default='0.1',
                      type="float",
                      help="Size of the plot point used to represent each pixel. [default: %default]")
    optionalGroup.add_option('-m','--map_res',
                      action="store",
                      dest="mapRes",
                      default='c',
                      type="choice",
                      choices=mapRes,
                      help="The map coastline resolution. Possible values are 'c' (coarse),'l' (low) and 'i' (intermediate). [default: %default]")
    optionalGroup.add_option('-a','--map_annotation',
                      action="store",
                      dest="mapAnn",
                      #default='',
                      type="string",
                      help="The map legend describing the dataset being shown. [default: IPPROD]")
    optionalGroup.add_option('-o','--output_file',
                      action="store",
                      dest="outputFile",
                      default="out.png",
                      type="string",
                      help="The full path of the output png file. [default: %default]")
    optionalGroup.add_option('--crefl_glob',
                      action="store",
                      dest="creflGlob",
                      type="string",
                      help="The full path of the VIIRS corrected relectance file output from crefl.")
    optionalGroup.add_option('--crefl_bin',
                      action="store",
                      dest="creflBin",
                      type="string",
                      help="The full path of the crefl binary executable.")
    optionalGroup.add_option('--crefl',
                      action="store_true",
                      dest="run_crefl",
                      help="Compute the corrected reflectances from VIIRS visible channels.")
    optionalGroup.add_option('--swath_image',
                      action="store_true",
                      dest="createSwathImage",
                      help="Produce an unprojected image file of the RGB scene.")
    optionalGroup.add_option('--no_map_image',
                      action="store_true",
                      dest="disableMapImage",
                      help="Do not produce the map-projected RGB image.")
    optionalGroup.add_option('--afires_file',
                      action="store",
                      dest="afilesFile",
                      type="string",
                      help="The full path of the VIIRS Active Fires file.")
    optionalGroup.add_option('-p','--product',
                      action="store",
                      dest="rgbProd",
                      type="choice",
                      default="TrueColor",
                      choices=prodChoices,
                      help='''The VIIRS RGB product.\n\n
                           Possible values are...
                           %s
                           ''' % (prodChoices.__str__()[1:-1]))
    optionalGroup.add_option('-v', '--verbose',
                      dest='verbosity',
                      action="count",
                      default=0,
                      help='each occurrence increases verbosity 1 level from ERROR: -v=WARNING -vv=INFO -vvv=DEBUG')

    parser.add_option_group(optionalGroup)

    # Parse the arguments from the command line
    (options, args) = parser.parse_args()

    # Set up the logging levels
    levels = [logging.ERROR, logging.WARN, logging.INFO, logging.DEBUG]
    output_logFile = "plot_viirs_RGB_crefl.py.log"

    try :
        logFormat = '(%(levelname)s): %(filename)s ; %(module)s::%(funcName)s ; line %(lineno)d -> %(message)s'
        #logging.basicConfig(format=logFormat, filename=output_logFile, filemode='w', level = levels[options.verbosity])
        logging.basicConfig(format=logFormat, level = levels[options.verbosity])

    except IndexError :
        print "ERROR : Invalid verbosity flag. Only -v, -vv, -vvv allowed, exiting..."
        sys.exit(1)

    verbosity = '-'
    for v in range(options.verbosity) :
        verbosity += 'v'

    # Check that the input files actually exist

    viirsGeo = False
    viirsRGB = False
    creflBin = False
    creflGlob = False

    if (options.geoFile is not None):
        if not path.exists(options.geoFile) :
            parser.error("Geolocation file \n\t%s\ndoes not exist, aborting..." % (options.geoFile))
        viirsGeo = True

    if (options.redFile is not None and \
        options.greenFile is not None and \
        options.blueFile is not None):
        if not path.exists(options.redFile) :
            parser.error("Channel file \n\t%s\ndoes not exist, aborting..." % (options.redFile))
        if not path.exists(options.greenFile) :
            parser.error("Channel file \n\t%s\ndoes not exist, aborting..." % (options.greenFile))
        if not path.exists(options.blueFile) :
            parser.error("Channel file \n\t%s\ndoes not exist, aborting..." % (options.blueFile))
        viirsRGB = True

    if (options.creflBin is not None):
        if not path.exists(str(options.creflBin)) :
            parser.error("crefl binary file \n\t%s\ndoes not exist, aborting..." % (options.creflBin))
        creflBin = True

    if (options.creflGlob is not None):
        creflList = glob(options.creflGlob)
        if (creflList == []) :
            parser.error("crefl file \n\t%s\ndoes not exist, aborting..." % (options.creflGlob))
        creflGlob = True

    # We have the VIIRS SDR files...
    if viirsRGB:

        LOG.info('VIIRS SDR files exist...')

        # Want to run crefl *and* have VIIRS geolocation...
        if options.run_crefl and creflBin and viirsGeo:

            if options.creflGlob is None :
                options.creflGlob = path.basename(options.geoFile)
                options.creflGlob = string.replace(string.replace(options.creflGlob,'GMTCO','CREFL'),'.h5','.hdf')
                options.creflGlob = path.join(path.curdir,options.creflGlob)
                LOG.info('Proposed output filename for crefl output: %s' % (options.creflGlob))
            
            LOG.info('Generating corrected reflectances using crefl...')
            _run_crefl(options.creflBin, options.geoFile, options.redFile, 
                    options.greenFile, options.blueFile,options.creflGlob) 

            LOG.info('Ingesting corrected reflectances generated using crefl...')
            RGBdict = gran_crefl(options.creflGlob,options.geoFile)

            LOG.info('Generating RGB map from corrected reflectances...')
            RGBdict = crefl_rgb(options.rgbProd,RGBdict,options.createSwathImage,swathName=options.outputFile)

        # Do *not* want to run crefl, *and* have VIIRS geolocation...
        elif not options.run_crefl and viirsGeo:

            LOG.info('Ingesting uncorrected reflectances...')
            #RGBdict = gran_Sdr(options.geoFile,options.redFile,options.greenFile,options.blueFile, \
                    #isRadiance=False)
            RGBdict = gran_Sdr2(options.geoFile,options.redFile,options.greenFile,options.blueFile, \
                    isRadiance=False)

            LOG.info('Generating RGB map from uncorrected reflectances...')
            LOG.info('Swathname = %s' %(options.outputFile))
            #RGBdict = viirs_rgb(options.rgbProd,RGBdict,options.createSwathImage)
            RGBdict = viirs_rgb(options.rgbProd,RGBdict, \
                    options.createSwathImage,swathName=options.outputFile)

    elif creflGlob :

        LOG.info('options.geoFile = %s' % (options.geoFile))
        creflList = glob(options.creflGlob)
        geoList = [] if options.geoFile is None else glob(options.geoFile)

        #LOG.info('creflList = %r' % (creflList))

        LOG.info('Ingesting corrected reflectances from existing crefl file %s' % (options.creflGlob))
        RGBdict = gran_crefl(creflList,geoList)

        sys.exit(0)

        LOG.info('Generating RGB map using corrected reflectances from existing crefl file %s' % (options.creflGlob))
        LOG.info('Swathname = %s' %(options.outputFile))
        RGBdict = crefl_rgb(options.rgbProd,RGBdict, \
                options.createSwathImage,swathName=options.outputFile)

    else :
        parser.error('''Incomplete mandatory arguments, must provide either
                red|green|blue files|geo|creflBin files, or
                geo|creflOut files, aborting...''')


    # Scale the three VIIRS channels, construct and RGB array
    #print "RGBdict['ModeGran'] = ",RGBdict['ModeGran']
    #if (RGBdict['ModeGran'] != 0):
        ##RGBdict = viirs_rgb(options.rgbProd,RGBdict)
        #RGBdict = crefl_rgb(options.rgbProd,RGBdict)
    #else :
        #pass

    #print "RGBdict.items() = ",RGBdict.items()

    # Make a plot if desired
    if not options.disableMapImage :

        prodFileName="" #'''%s\n%s REV:%s''' % (path.basename(options.geoFile),options.svnRepoPath,str(options.svnRevision))
        dataSet = string.lower((options.rgbProd).lstrip())
        mapRes = str(options.mapRes)

        # Copy the Channel name to the map annotation
        if options.mapAnn == None :
            options.mapAnn = options.rgbProd
        else :
            pass
            options.mapAnn = string.replace(options.mapAnn,'\\n','\n')

        _create_map_image_RGB(RGBdict,pointSize=options.pointSize,scale=options.scale,mapRes=mapRes,stride=options.stride,\
                sdrKey='M15',isRadiance=False,prodFileName=prodFileName,\
                outFileName=options.outputFile,dpi=options.dpi,titleStr=options.mapAnn)
        #orthoPlot_RGB_noravel(RGBdict,pointSize=options.pointSize,scale=options.scale,mapRes=mapRes,stride=options.stride,\
                #sdrKey='M15',isRadiance=False,prodFileName=prodFileName,\
                #outFileName=options.outputFile,dpi=options.dpi,titleStr=options.mapAnn)


    print "Exiting..."
    sys.exit(0)

if __name__ == '__main__':
    main()