-
Geoff Cureton authoredgit-svn-id: https://svn.ssec.wisc.edu/repos/geoffc/Python@311 f45efd6a-04c6-4179-886d-a06c9f63a260
Geoff Cureton authoredgit-svn-id: https://svn.ssec.wisc.edu/repos/geoffc/Python@311 f45efd6a-04c6-4179-886d-a06c9f63a260
plot_viirs_RGB_crefl.py 59.21 KiB
#!/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()