diff --git a/pyglance/glance/imapp_plot.py b/pyglance/glance/imapp_plot.py
index 4db9aa4d0a11f1f35b6b74ece76bf065cefe4dce..bc40328f7fecdb3df37fb01992369246a5e208ea 100644
--- a/pyglance/glance/imapp_plot.py
+++ b/pyglance/glance/imapp_plot.py
@@ -20,6 +20,7 @@ import matplotlib.pyplot as plt
import matplotlib.colors as colors
from mpl_toolkits.basemap import Basemap
+import re
import numpy as np
import glance.data as dataobj
@@ -27,19 +28,31 @@ import glance.data as dataobj
LOG = logging.getLogger(__name__)
defaultValues = {
- 'longitudeVar': 'xtraj',
- 'latitudeVar': 'ytraj',
- 'initAODVar': 'aod_traj',
- 'trajPressVar': 'ptraj',
- 'timeVar': 'time',
- 'nePiece': 'NE_',
- 'swPiece': 'SW_',
- 'latPiece': 'LAT',
- 'lonPiece': 'LON',
- 'figureName': 'frame.png',
- 'figureDPI': 200
+ 'longitudeVar': 'xtraj',
+ 'latitudeVar': 'ytraj',
+ 'initAODVar': 'aod_traj',
+ 'trajPressVar': 'ptraj',
+ 'timeVar': 'time',
+ 'neName': 'TOP_RIGHT_LON_LAT', #TODO need to fix this? 'NE',
+ 'swName': 'BOTTOM_LEFT_LON_LAT', # TODO need to fix this? 'SW',
+ 'windLonName': 'Longitude',
+ 'windLatName': 'Latitude',
+ 'windTimeName': 'time',
+ 'windUName': 'uwind',
+ 'windVName': 'vwind',
+ 'emissNameBase': 'Cloud_Effective_Emissivity_',
+ 'aodNameBase': 'Optical_Depth_Land_And_Ocean_',
+ 'optLonBase': 'Lon_',
+ 'optLatBase': 'Lat_',
+ 'figureName': 'frame.png',
+ 'figureDPI': 200
}
+ACCEPTABLE_MAP_PROJECTIONS = ['cylindrical']
+
+LEVELS_FOR_BACKGROUND_PLOTS = 20
+DEFAULT_WINDS_SUBDIVIDE_FACTOR = 1
+
# a custom colormap or the Trajectory Pressures
color_data = {
'red' : ( (0.0, 1.0, 1.0),
@@ -72,70 +85,299 @@ color_data = {
}
light_trajectory_pressure_color_map = matplotlib.colors.LinearSegmentedColormap('lightTrajPressCM', color_data, 256)
-def _create_imapp_figure (initAODdata, initLongitudeData, initLatitudeData,
- pressureData=None, pressLongitudeData=None, pressLatitudeData=None,
- baseMapInstance=None, figureTitle="MODIS AOD & AOD Trajectories",
- useDarkBackground=True, latRange=None, lonRange=None) :
+def _check_requested_projection (projection, doWarn=True, fileDescription="") :
"""
- TODO, I'm still deciding how this funciton works so it will probably change drastically
+ Check the requested projection, issuing a warning if it is not available and doWarn is true
"""
- # build the plot
- figure = plt.figure()
- tempColor = 'k' if useDarkBackground else 'w' # either use a black or white background
- axes = figure.add_subplot(111, axisbg=tempColor)
+ LOG.debug("Requested map projection: " + str(projection))
+ if projection not in ACCEPTABLE_MAP_PROJECTIONS :
+ LOG.warn("The requested map projection in " + fileDescription + " (" + str(projection)
+ + ") is not a projection that can be produced by this program. "
+ + "A default projection will be produced. Output may not be displayed as expected.")
+
+def _get_matching_terms_from_list (list_of_strings, regular_expression_text, case_sensitive=False) :
+ """
+ given a list and a text regular expression in the form taken by the python re module,
+ return a list of all strings in the list that match that regular expression
- # build extra info to go to the map plotting function
- kwargs = { }
+ TODO, there's probably a better way to do this, but this is what I've got for now
+ """
- # draw the basic physical and geopolitical features
- tempColor='w' if useDarkBackground else 'k' # either draw black or white lines
- if baseMapInstance is not None :
- baseMapInstance.drawcoastlines( color=tempColor, linewidth=0.5)
- baseMapInstance.drawcountries( color=tempColor, linewidth=0.5)
- baseMapInstance.drawstates( color=tempColor, linewidth=0.5)
- baseMapInstance.drawmapboundary(color=tempColor, linewidth=0.5)
- # draw the parallels and meridians
- if latRange is not None :
- parallels = arange(-80., 90., latRange / 4.0)
- baseMapInstance.drawparallels(parallels,labels=[1,0,0,1], color=tempColor, linewidth=0.5)
- if lonRange is not None :
- meridians = arange(0., 360., lonRange / 4.0)
- baseMapInstance.drawmeridians(meridians,labels=[1,0,0,1], color=tempColor, linewidth=0.5)
+ tempExpression = re.compile(regular_expression_text) if case_sensitive else re.compile(regular_expression_text, re.IGNORECASE)
- # translate the longitude and latitude data sets into map coordinates
- pressX, pressY = None, None
- if pressureData is not None :
- pressX, pressY = baseMapInstance(pressLongitudeData, pressLatitudeData)
- initX, initY = baseMapInstance(initLongitudeData, initLatitudeData)
+ toReturn = [ ]
+ for term in list_of_strings :
+ if tempExpression.match(term) is not None :
+ toReturn.append(term)
+ return toReturn
+
+def _find_most_current_index_at_time (times_array, current_time) :
+ """
+ figure out the index of the data in times_array that is most
+ current given that current_time is now
+
+ if for some reason the data in the times_array is all after
+ current_time, return -1
+ otherwise, return the index of the most current data that
+ has happened or is happening now (no future indexes!)
+
+ note: it is assumed that the times_array is arranged in
+ increasing order (earlier to later times) and no time is
+ repeated.
+ """
+
+ toReturn = -1
+
+ for indexNum in range(len(times_array)) :
+ tempTime = times_array[indexNum]
+ if tempTime <= current_time :
+ toReturn = indexNum
+
+ return toReturn
+
+def _draw_contour_with_basemap (baseMapInstance, data, lonData, latData, levels=None, **kwargs) :
+ """
+ draw a contour plot of the data using the basemap and the provided lon and lat
+ """
+
+ # translate into the coordinate system of the basemap
+ tempX, tempY = baseMapInstance(lonData, latData)
+
+ # show the plot if there is data
+ if (data is not None) :
+
+ if levels is not None :
+ p = baseMapInstance.contourf(tempX, tempY, data, levels, **kwargs)
+ else :
+ p = baseMapInstance.contourf(tempX, tempY, data, **kwargs)
+
+ # TODO, do I need to return tempX and tempY in the future?
+
+def _draw_winds_with_basemap (baseMapInstance,
+ uData, vData,
+ lonData, latData,
+ valueData=None,
+ **kwargs) :
+ """
+ draw the given wind vectors using the provided basemap
+
+ if the value data isn't None, use it to control the color values plotted on the vectors
+ """
+
+ # translate into the coordinate system of the basemap
+ tempX, tempY = baseMapInstance(lonData, latData)
+
+ # show the quiver plot if there is data
+ if (uData is not None) and (vData is not None) :
+
+ # TODO, it looks like there's a bug that can affect the data sent into quiver, for now copy it
+ # bug report can be found here: https://github.com/matplotlib/matplotlib/issues/625
+ tempUData = uData.copy()
+ tempVData = vData.copy()
+
+ kwargs['width'] = 0.001
+ kwargs['headwidth'] = 4
+
+ if valueData is None:
+ p = baseMapInstance.quiver(tempX, tempY, tempUData, tempVData, **kwargs)
+ else :
+ p = baseMapInstance.quiver(tempX, tempY, tempUData, tempVData, valueData, **kwargs)
+
+ # TODO, do I need to return tempX and tempY in the future?
+
+def _plot_pressure_data (baseMapInstance, pressureLat, pressureLon, pressureData=None, colorMap=dark_trajectory_pressure_color_map) :
+ """
+ plot the pressure data and return the color bar so it can be moved around as needed
+ """
+
+ colorBarToReturn = None
- color_map_to_use = dark_trajectory_pressure_color_map if useDarkBackground else light_trajectory_pressure_color_map
# if we have pressure data plot that
if pressureData is not None :
+
+ # translate the longitude and latitude into map coordinates
+ pressX, pressY = baseMapInstance(pressureLon, pressureLat)
+
# I'm taking advantage of the fact that I can remove the black edge line with lw=0 (line width = 0) and scale the marker down with s=0.5
- baseMapInstance.scatter(pressX, pressY, s=0.5, c=pressureData, marker='o', cmap=color_map_to_use, vmin=0, vmax=1000, lw=0)
+ baseMapInstance.scatter(pressX, pressY, s=3.0, c=pressureData, marker='o', cmap=colorMap, vmin=0, vmax=1000, lw=0)
# create the pressure colorbar
- cbar1 = colorbar(format='%.3g', orientation='horizontal', shrink=0.25)
- cbar1.ax.set_position([0.1, -0.16, 0.25, 0.25])
- for tempText in cbar1.ax.get_xticklabels():
- #print(tempText)
+ colorBarToReturn = colorbar(format='%.3g', orientation='horizontal', shrink=0.25)
+ for tempText in colorBarToReturn.ax.get_xticklabels():
tempText.set_fontsize(5)
- cbar1.set_label("Trajectory Pressure (mb)")
+ colorBarToReturn.set_label("Trajectory Pressure (mb)")
+
+ return colorBarToReturn
+
+def _plot_initial_modis_aod (baseMapInstance, longitude, latitude, initAODdata, colorMap=cm.jet) :
+ """
+ plot initial modis aod points using the provided base map, return the created
+ color bar so it can be moved around as needed
+ """
- # plot the origin points after the pressure so they'll appear on top, I'm assuming we will always have this data
- baseMapInstance.scatter(initX, initY, s=10, c=initAODdata, marker='o', cmap=cm.jet, vmin=0.0, vmax=1.0, lw=0.5)
+ # translate the longitude and latitude into map coordinates
+ initX, initY = baseMapInstance(longitude, latitude)
+
+ # plot the origin points, these are being scaled to have a thin but visible black border line
+ baseMapInstance.scatter(initX, initY, s=10, c=initAODdata, marker='o', cmap=colorMap, vmin=0.0, vmax=1.0, lw=0.5)
# make a color bar
- cbar2 = colorbar(format='%.3g', orientation='horizontal', shrink=0.25)
- cbar2.ax.set_position([0.4, -0.16, 0.25, 0.25])
- for tempText in cbar2.ax.get_xticklabels():
- #print(tempText)
+ colorBarToReturn = colorbar(format='%.3g', orientation='horizontal', shrink=0.25)
+ for tempText in colorBarToReturn.ax.get_xticklabels():
tempText.set_fontsize(5)
- cbar2.set_label("MODIS AOD") # TODO, how do I get a second Trajectory Pressure colorbar in the right place?
+ colorBarToReturn.set_label("MODIS AOD")
+
+ return colorBarToReturn
+
+def _build_basic_figure_with_map (baseMapInstance, latRange=None, lonRange=None, useDarkBackground=True,) :
+ """
+ create a figure with the longitude and latitude info given
+ """
+
+ # create the empty figure
+ figure = plt.figure()
+ bkgdColor = 'k' if useDarkBackground else 'w' # either use a black or white background
+ axes = figure.add_subplot(111, axisbg=bkgdColor)
+
+ # draw the basic physical and geopolitical features
+ lineColor = 'w' if useDarkBackground else 'k' # either draw black or white lines
+ if baseMapInstance is not None :
+ baseMapInstance.drawcoastlines( color=lineColor, linewidth=0.5)
+ baseMapInstance.drawcountries( color=lineColor, linewidth=0.5)
+ baseMapInstance.drawstates( color=lineColor, linewidth=0.5)
+ baseMapInstance.drawmapboundary(color=lineColor, linewidth=0.5)
+ # draw the parallels and meridians
+ if latRange is not None :
+ parallels = arange(-80., 90., latRange / 4.0)
+ baseMapInstance.drawparallels(parallels,labels=[1,0,0,1], color=lineColor, linewidth=0.5)
+ if lonRange is not None :
+ meridians = arange( 0., 360., lonRange / 4.0)
+ baseMapInstance.drawmeridians(meridians,labels=[1,0,0,1], color=lineColor, linewidth=0.5)
+
+ return axes, figure
+
+# todo, I hate this solution for making contourf show the number of levels that I want,
+# but I've searched repeatedly and can't find another way to get it to let me specify
+def _make_range_with_data_objects(list_of_data_objects, num_intervals, offset_to_range=0.0) :
+ """
+ get an array with numbers representing the bounds of a set of ranges
+ that covers all the valid data present in the data objects given
+ (these may be used for plotting the data)
+ if an offset is passed, the outtermost range will be expanded by that much
+
+ note: the list of data objects may be 0 or more data objects, if it 0, this
+ method returns None since no useful range can be created
+ """
+
+ if len(list_of_data_objects) <= 0 :
+ return None
+
+ minVal = np.PINF # positive infinity
+ maxVal = np.NINF # negative infinity
+
+ # include each data set in the min/max calculations
+ for each_data_object in list_of_data_objects :
+ minVal = min(each_data_object.get_min(), minVal)
+ maxVal = max(each_data_object.get_max(), maxVal)
+
+ # TODO, there's the possibility for failure here if all data is fully masked out?
+
+ # take any offsets into account
+ minVal = minVal - offset_to_range
+ maxVal = maxVal + offset_to_range
+
+ return np.linspace(minVal, maxVal, num_intervals)
+
+# todo, I hate this solution for making contourf show the number of levels that I want,
+# but I've searched repeatedly and can't find another way to get it to let me specify
+def _make_range_with_masked_arrays(list_of_masked_arrays, num_intervals, offset_to_range=0.0) :
+ """
+ get an array with numbers representing the bounds of a set of ranges
+ that covers all the valid data present in the masked arrays given
+ (these may be used for plotting the data)
+ if an offset is passed, the outtermost range will be expanded by that much
+
+ note: the list of masked arrays may be 0 or more masked arrays, if it 0, this
+ method returns None since no useful range can be created
+ """
+
+ if len(list_of_masked_arrays) <= 0 :
+ return None
+
+ minVal = np.PINF # positive infinity
+ maxVal = np.NINF # negative infinity
+
+ # include each data set in the min/max calculations
+ for each_masked_array in list_of_masked_arrays :
+ minVal = min(each_masked_array.min(), minVal)
+ maxVal = max(each_masked_array.max(), maxVal)
+
+ # TODO, there's the possibility for failure here if all data is fully masked out?
+
+ # take any offsets into account
+ minVal = minVal - offset_to_range
+ maxVal = maxVal + offset_to_range
+
+ return np.linspace(minVal, maxVal, num_intervals)
+
+def _create_imapp_figure (initAODdata, initLongitudeData, initLatitudeData,
+ pressureData=None, pressLongitudeData=None, pressLatitudeData=None,
+ baseMapInstance=None, figureTitle="MODIS AOD & AOD Trajectories",
+ useDarkBackground=True, latRange=None, lonRange=None,
+ windsDataU=None, windsDataV=None, windsDataLon=None, windsDataLat=None,
+ backgroundDataSets={ }) :
+ """
+ TODO, I'm still deciding how this funciton works so it may change
+
+ backgroundDataSets should be a dictionary in the form
+ {
+ 1: [dataToContourPlot1, LatitudeData1, LongitudeData1, colorMap1, levels],
+ 2: [dataToContourPlot2, LatitudeData2, LongitudeData2, colorMap2, levels],
+ ...
+ N: [dataToContourPlotN, LatitudeDataN, LongitudeDataN, colorMapN],
+ }
+ The keys should all be numbers and the data sets will be plotted in order from lowest key value to highest (allowing z-order control).
+
+ """
+
+ # create a figure and draw geopolitical features on it
+ axes, figure = _build_basic_figure_with_map (baseMapInstance,
+ latRange=latRange, lonRange=lonRange,
+ useDarkBackground=useDarkBackground)
+
+ # choose the color map
+ color_map_to_use = dark_trajectory_pressure_color_map if useDarkBackground else light_trajectory_pressure_color_map
+
+ # plot out any background data sets that we have
+ for orderKey in sorted(backgroundDataSets.keys()) :
+ tempDataSet, tempLatitude, tempLongitude, tempColorMap, tempLevelsList = backgroundDataSets[orderKey]
+ _draw_contour_with_basemap (baseMapInstance, tempDataSet, tempLongitude, tempLatitude, cmap=tempColorMap, levels=tempLevelsList, lw=0, antialiased=False)
+ # lw is line width, use 0 to hide the lines between contours
+ # there appears to be a bug in sub-pixel anti-aliasing that is making
+ # contour lines show up even with lw=0, so turn that off for now :(
+ # see also: http://stackoverflow.com/questions/8263769/hide-contour-linestroke-on-pyplot-contourf-to-get-only-fills
+
+ # if we have winds to draw, draw those first so they're on the bottom under all the other data
+ if ((windsDataU is not None) and (windsDataV is not None) and (windsDataLon is not None) and (windsDataLat is not None)) :
+ tempColor = 'w' if useDarkBackground else 'k' # plot our winds in a color that will show up on the background
+ _draw_winds_with_basemap (baseMapInstance, windsDataU, windsDataV, windsDataLon, windsDataLat, color=tempColor)
+
+ # plot the pressure data if appropriate
+ pressColorBar = _plot_pressure_data(baseMapInstance, pressLatitudeData, pressLongitudeData, pressureData=pressureData, colorMap=color_map_to_use)
+ # if we got a color bar back, make sure it's in the right place
+ if pressColorBar is not None :
+ pressColorBar.ax.set_position([0.4, -0.16, 0.25, 0.25])
+
+ # plot the initial modis AOD points after the pressure so the AOD points are always on top (and visible)
+ aodColorBar = _plot_initial_modis_aod(baseMapInstance, initLongitudeData, initLatitudeData, initAODdata)
+ # if we got a color bar back, make sure it's in the right place
+ if aodColorBar is not None :
+ aodColorBar.ax.set_position([0.1, -0.16, 0.25, 0.25])
# now that we've moved everything around, make sure our main image is in the right place
+ # TODO, make sure this resizing/placement will work in more general cases
axes.set_position([0.1, 0.15, 0.8, 0.8]) # why was this method so hard to find?
# set up the figure title
@@ -151,7 +393,7 @@ def main():
run "%prog help" to list commands
examples:
-python -m glance.imapp_plot plot A.nc
+python -m glance.imapp_plot aodTraj A.nc
"""
# the following represent options available to the user on the command line:
@@ -166,30 +408,38 @@ python -m glance.imapp_plot plot A.nc
parser.add_option('-w', '--debug', dest="debug",
action="store_true", default=False, help="enable debug output")
- # file generation related options TODO, make this work
+ # file generation related options
parser.add_option('-p', '--outpath', dest="outpath", type='string', default='./',
help="set path to output directory")
- # file variable settings
+ # file variable settings TODO, do we need these?
parser.add_option('-o', '--longitude', dest="longitudeVar", type='string',
help="set name of longitude variable")
parser.add_option('-a', '--latitude', dest="latitudeVar", type='string',
help="set name of latitude variable")
+ # display related settings
+ parser.add_option('-d', '--windsSubdivideFactor', dest="subdivideFactor", type='int',
+ default=1, help="factor to subdivide and filter out some winds data for visiblity; " +
+ "higher factors result in fewer winds; data is filtered evenly along the indices of the winds data grid")
+
# time related settings
parser.add_option('-s', '--startTime', dest="startTime", type='int',
default=0, help="set first time to process")
parser.add_option('-e', '--endTime', dest="endTime", type='int',
help="set last time to process")
parser.add_option('-f', '--futureWindow', dest="futureWindow", type='int',
- default=6, help="set number of hours of future pressures to show")
+ default=6, help="set number of hours of future pressures to show; 6 hours is the default")
+ parser.add_option('-j', '--jumpEmptyStartHours', dest="doJump",
+ action="store_true", help="start generating images on the first hour with trajectory data that is after 0 (ie. not 0); if used this flag " +
+ "overrrides the start time option when determining which images will be created at the start of a series")
parser.add_option('-t', '--test', dest="self_test",
action="store_true", default=False, help="run internal unit tests")
# TODO will add this once we're out of alpha
#parser.add_option('-n', '--version', dest='version',
- # action="store_true", default=False, help="view the glance version")
+ # action="store_true", default=False, help="view the imapp plot version")
# parse the uers options from the command line
@@ -218,13 +468,18 @@ python -m glance.imapp_plot plot A.nc
The following functions represent available menu selections.
"""
- def plot(*args):
- """plot trajectory frames
+ def aodTraj(*args):
+ """plot AOD trajectory frames
Given a file with trajectory possitions and pressures over time, plot out
images of these trajectories on the Earth and save it to disk.
+
+ optionally, if a second file with winds data and MODIS cloud effective
+ emissivity data is provided, this will also be plotted at comprable time
+ steps for each frame
"""
LOG.debug("startTime: " + str(options.startTime))
+ LOG.debug("will startTime be modified by jumping frames? " + str(options.doJump))
LOG.debug("endTime: " + str(options.endTime))
LOG.debug("futureWindow: " + str(options.futureWindow))
@@ -253,23 +508,126 @@ python -m glance.imapp_plot plot A.nc
trajectoryTimeData = trajectoryFileObject.file_object[defaultValues['timeVar']]
# get information on where we should display the data
- northeastLon = trajectoryFileObject.file_object.get_global_attribute( defaultValues['nePiece'] + defaultValues['lonPiece'] )
- northeastLat = trajectoryFileObject.file_object.get_global_attribute( defaultValues['nePiece'] + defaultValues['latPiece'] )
- southwestLon = trajectoryFileObject.file_object.get_global_attribute( defaultValues['swPiece'] + defaultValues['lonPiece'] )
- southwestLat = trajectoryFileObject.file_object.get_global_attribute( defaultValues['swPiece'] + defaultValues['latPiece'] )
- latRange = abs(northeastLat - southwestLat)
- lonRange = abs(southwestLon - northeastLon)
+ northeastLon, northeastLat = trajectoryFileObject.file_object.get_global_attribute( defaultValues['neName'] )
+ southwestLon, southwestLat = trajectoryFileObject.file_object.get_global_attribute( defaultValues['swName'] )
+ latRange = abs(northeastLat - southwestLat)
+ lonRange = abs(southwestLon - northeastLon)
+ LOG.debug ("Viewing window, northeast corner (lat / lon): " + str(northeastLat) + " / " + str(northeastLon))
+ LOG.debug ("Viewing window, southwest corner (lat / lon): " + str(southwestLat) + " / " + str(southwestLon))
+
+ # double check the map projection
+ _check_requested_projection(trajectoryFileObject.file_object.get_global_attribute("MAP_PROJECTION"), fileDescription="the main trajectory file")
+
+ # check to see if we have a second file for the optional data
+ LOG.info("Attempting to open optional file (if present)")
+ doOptionalWinds, doOptionalCloudEffectiveEmiss, doOptionalDepthLandAndOcean = False, False, False
+ optionalDataWindsLatitude, optionalDataWindsLongitude, optionalDataWindsTime, optionalDataWindU, optionalDataWindV = None, None, None, None, None
+ listOfCloudEffectiveEmissivityPasses, listOfOpticalDepthLandAndOcean, listOfOptionalLon, listOfOptionalLat = { }, { }, { }, { }
+ if len(args) > 1 :
+ optionalFilePath = str(args[1])
+ LOG.debug("Opening optional file: " + optionalFilePath)
+ optionalFileObject = dataobj.FileInfo(optionalFilePath)
+ if optionalFileObject is None :
+ LOG.warn("Optional file (" + optionalFilePath + ") could not be opened.")
+ LOG.warn("Optional winds and cloud effective emissivity data will not be loaded/plotted.")
+ else :
+
+ # hang on to the list of variables for some checks later
+ tempVarList = optionalFileObject.file_object()
+
+ # get info on where this data wants to be displayed
+ optionalNELon, optionalNELat = optionalFileObject.file_object.get_global_attribute( defaultValues['neName'] )
+ optionalSWLon, optionalSWLat = optionalFileObject.file_object.get_global_attribute( defaultValues['swName'] )
+ LOG.debug ("Viewing window in optional file, northeast corner (lat / lon): " + str(optionalNELat) + " / " + str(optionalNELon))
+ LOG.debug ("Viewing window in optional file, southwest corner (lat / lon): " + str(optionalSWLat) + " / " + str(optionalSWLon))
+
+ # if the two viewing windows are not the same, we don't want to display this data TODO, confirm that this is the right strategy?
+ if ((optionalNELat != northeastLat) or (optionalNELon != northeastLon) or
+ (optionalSWLat != southwestLat) or (optionalSWLon != southwestLon)) :
+ LOG.debug ("Incompatable viewing area given in optional file. Optional data will not be displayed.")
+
+ else :
+
+ # double check the map projection
+ _check_requested_projection(optionalFileObject.file_object.get_global_attribute("MAP_PROJECTION"), fileDescription="the optional data file")
+
+ # otherwise try to get the data we need if possible
+ if ((defaultValues['windLonName'] in tempVarList) and (defaultValues['windLatName'] in tempVarList) and
+ (defaultValues['windTimeName'] in tempVarList) and
+ (defaultValues['windUName'] in tempVarList) and (defaultValues['windVName'] in tempVarList)) :
+
+ # if we have all the winds related variables available, load the winds data
+ # TODO, should the user be able to control these names?
+ optionalDataWindsLatitude = optionalFileObject.file_object[defaultValues['windLatName']]
+ optionalDataWindsLongitude = optionalFileObject.file_object[defaultValues['windLonName']]
+ optionalDataWindsTime = optionalFileObject.file_object[defaultValues['windTimeName']]
+ optionalDataWindU = optionalFileObject.file_object[defaultValues['windUName']]
+ optionalDataWindV = optionalFileObject.file_object[defaultValues['windVName']]
+ doOptionalWinds = True
+
+ # TODO, I'd really rather not do this, for the moment I don't have a choice
+ tempLatSize = optionalDataWindsLatitude.size
+ tempLonSize = optionalDataWindsLongitude.size
+ tempLat = optionalDataWindsLatitude
+ tempLon = optionalDataWindsLongitude
+ optionalDataWindsLatitude = np.transpose(np.repeat([tempLat], tempLonSize, axis=0))
+ optionalDataWindsLongitude = np.repeat([tempLon], tempLatSize, axis=0)
+
+ # if any cloud effective emissivity swaths are available load them
+ possibleEmissNames = _get_matching_terms_from_list(tempVarList, defaultValues['emissNameBase'])
+ if len(possibleEmissNames) > 0 :
+ LOG.debug("cloud effective emissivity variables found: " + str(possibleEmissNames))
+ for emissName in possibleEmissNames :
+ tempEmissData = optionalFileObject.file_object[emissName]
+ tempMissingValue = optionalFileObject.file_object.missing_value(emissName)
+ listOfCloudEffectiveEmissivityPasses[emissName] = ma.masked_where(tempEmissData == tempMissingValue, tempEmissData)
+ doOptionalCloudEffectiveEmiss = True
+
+ # load optical depth land and ocean if it's present
+ possibleOpticalDepthNames = _get_matching_terms_from_list(tempVarList, defaultValues['aodNameBase'])
+ if len(possibleOpticalDepthNames) > 0 :
+ LOG.debug("optical depth land and ocean variables found: " + str(possibleOpticalDepthNames))
+ for aodName in possibleOpticalDepthNames :
+ tempAODData = optionalFileObject.file_object[aodName]
+ tempMissingValue = optionalFileObject.file_object.missing_value(aodName)
+ listOfOpticalDepthLandAndOcean[aodName] = ma.masked_where(tempAODData == tempMissingValue, tempAODData)
+ doOptionalDepthLandAndOcean = True
+
+ # if either the cloud effective emissivity or the optical depth land and ocean were loaded,
+ # then we need to load the associated lat and lon variables
+ if doOptionalCloudEffectiveEmiss or doOptionalDepthLandAndOcean :
+ possibleOptionalLatNames = _get_matching_terms_from_list(tempVarList, defaultValues['optLatBase'])
+ possibleOptionalLonNames = _get_matching_terms_from_list(tempVarList, defaultValues['optLonBase'])
+ LOG.debug ("optional longitude variables found: " + str(possibleOptionalLonNames))
+ LOG.debug ("optional latitude variables found: " + str(possibleOptionalLatNames))
+ # if we can't find lon/lat, we can't plot these! Note: this still fails if we have partial lon/lat, could compare sizes in FUTURE
+ if (len(possibleOptionalLatNames) <= 0) or (len(possibleOptionalLonNames) <= 0) :
+ LOG.warn("Unable to find corresponding longitude and latitude data for some optional variables. "
+ + "Data without corresponding longitude and latitude data will not be plotted.")
+ doOptionalDepthLandAndOcean = False
+ doOptionalCloudEffectiveEmiss = False
+ else :
+ # if we found lon and lat, load them
+ for lonName in possibleOptionalLonNames :
+ listOfOptionalLon[lonName] = optionalFileObject.file_object[lonName]
+ for latName in possibleOptionalLatNames :
+ listOfOptionalLat[latName] = optionalFileObject.file_object[latName]
# build a basemap
LOG.info("Building basemap object.")
- projectionName = 'merc' # use the Mercator Projection
+ projectionName = 'merc' # use the Mercator Projection; TODO at some point this will need to be checked with a global attribute
basemapObject = Basemap (projection=projectionName,llcrnrlat=southwestLat,urcrnrlat=northeastLat,
llcrnrlon=southwestLon, urcrnrlon=northeastLon, lat_ts=20, resolution='l') # TODO do I need to use lat_ts=20, ?
# sort out the times we're using
futureWindow = options.futureWindow
- startTime = options.startTime
- endTime = options.endTime if options.endTime is not None else trajectoryTimeData[-1]
+ startTime = options.startTime if options.startTime >= 0 else 0
+ # if we need to jump to start at the first trajectory data after 0, do that now
+ if options.doJump :
+ startTimeTemp = trajectoryTimeData[0] if trajectoryTimeData[0] > 0 else trajectoryTimeData[1]
+ startTime = max(startTime, startTimeTemp)
+ # I am assuming that the time data will never be negative
+ endTime = options.endTime if options.endTime is not None else trajectoryTimeData[-1]
# as a sanity check, it's not really productive to make frames for times after our data ends
if endTime > trajectoryTimeData[-1] :
endTime = trajectoryTimeData[-1]
@@ -280,41 +638,95 @@ python -m glance.imapp_plot plot A.nc
initLatData = latitudeData[0].ravel()
for currentTime in range (startTime, endTime + 1) :
- # select only the window of trajectory data we need
- alreadyHappenedTimeIndex = 0
- tooFarInFutureTimeIndex = trajectoryTimeData.size
- for tempIndex in range (0, trajectoryTimeData.size) :
-
- # first find the time that corresponds to the "most recently happened" index
- if trajectoryTimeData[tempIndex] <= currentTime :
- # TODO if time doesn't always increase, this needs another check
- alreadyHappenedTimeIndex = tempIndex
-
- # then figure out how much further we can go before we pass the futureWindow's edge
- if (trajectoryTimeData[tempIndex] > (currentTime + futureWindow)) :
- # TODO, if time data doesn't always increase I also need & (trajectoryTimeData[tempIndex] < trajectoryTimeData[tooFarInFutureTimeIndex])
- tooFarInFutureTimeIndex = tempIndex
- break # TODO, this assumes time data always increases; also breaks suck so eventually take this out
+ # TEMP for debugging during develpment
+ LOG.debug ("**********")
+
+ # get the time that represents the most current trajectory data
+ mostCurrentTrajTimeIndex = _find_most_current_index_at_time (trajectoryTimeData, currentTime)
+ # get the time that represents the last point included in the future window
+ lastTimeIncludedIndex = _find_most_current_index_at_time (trajectoryTimeData, (currentTime + futureWindow))
# only get data to plot the trajectory points if there is some available
thisFramePressures = None
thisFramePressLon = None
thisFramePressLat = None
- if alreadyHappenedTimeIndex + 1 < tooFarInFutureTimeIndex :
+ if mostCurrentTrajTimeIndex <= lastTimeIncludedIndex :
- LOG.debug("Already happened index: " + str(alreadyHappenedTimeIndex))
- LOG.debug("Too far future index: " + str(tooFarInFutureTimeIndex))
+ # TODO , double check that this isn't off by one
+ LOG.debug("Most current time index: " + str(mostCurrentTrajTimeIndex))
+ LOG.debug("Last data in future window index: " + str(lastTimeIncludedIndex))
+ LOG.debug("Most current time: " + str(trajectoryTimeData[mostCurrentTrajTimeIndex]))
+ LOG.debug("Last data in future window time: " + str(trajectoryTimeData[lastTimeIncludedIndex]))
# now pull out the pressure data and related lon/lat info for plotting
- thisFramePressures = trajectoryPressureData[alreadyHappenedTimeIndex + 1 : tooFarInFutureTimeIndex].ravel()
- thisFramePressLon = longitudeData[alreadyHappenedTimeIndex + 1 : tooFarInFutureTimeIndex].ravel()
- thisFramePressLat = latitudeData[alreadyHappenedTimeIndex + 1 : tooFarInFutureTimeIndex].ravel()
-
+ thisFramePressures = trajectoryPressureData[mostCurrentTrajTimeIndex : lastTimeIncludedIndex+1].ravel()
+ thisFramePressLon = longitudeData[mostCurrentTrajTimeIndex : lastTimeIncludedIndex+1].ravel()
+ thisFramePressLat = latitudeData[mostCurrentTrajTimeIndex : lastTimeIncludedIndex+1].ravel()
+
+ # we need to plot the winds under the other data so they don't cover it up
+ currentWindsU, currentWindsV, currentWindsLon, currentWindsLat = None, None, None, None
+ if doOptionalWinds :
+ # figure out which winds data we should use
+ mostCurrentWindsTimeIndex = _find_most_current_index_at_time (optionalDataWindsTime, currentTime)
+ LOG.debug("Optional winds most current time index: " + str(mostCurrentWindsTimeIndex))
+ LOG.debug("Time of most recent optional winds: " + str(optionalDataWindsTime[mostCurrentWindsTimeIndex]))
+ # check to make sure we found any winds at the current time
+ if mostCurrentWindsTimeIndex > -1 :
+ currentWindsU = optionalDataWindU[mostCurrentWindsTimeIndex]
+ currentWindsV = optionalDataWindV[mostCurrentWindsTimeIndex]
+
+ maskingFactor = options.subdivideFactor if options.subdivideFactor > 0 else DEFAULT_WINDS_SUBDIVIDE_FACTOR
+ LOG.debug ("Using winds subdivide factor: " + str(maskingFactor))
+ maskTemp = np.zeros(currentWindsV.shape, dtype=np.bool)
+ if maskingFactor > 1 :
+ maskTemp1 = np.ones(currentWindsV.shape, dtype=np.bool)
+ maskTemp2 = np.ones(currentWindsV.shape, dtype=np.bool)
+ numbers1 = np.repeat([range(currentWindsV.shape[1])], currentWindsV.shape[0], axis=0)
+ numbers2 = np.transpose(np.repeat([range(currentWindsV.shape[0])], currentWindsV.shape[1], axis=0))
+ maskTemp1 [ numbers1 % maskingFactor == 0 ] = False
+ maskTemp2 [ numbers2 % maskingFactor == 0 ] = False
+ maskTemp = maskTemp1 | maskTemp2
+ currentWindsU = np.ma.array(currentWindsU, mask=maskTemp)
+ currentWindsV = np.ma.array(currentWindsV, mask=maskTemp)
+
+ # build up the background data sets we want to plot
+ backgroundData = { }
+ zCounter = 1
+
+ # if we're doing cloud effective emissivity, add those to the background list
+ if doOptionalCloudEffectiveEmiss :
+ tempLevels = _make_range_with_masked_arrays(listOfCloudEffectiveEmissivityPasses.values(), LEVELS_FOR_BACKGROUND_PLOTS)
+ for emissName in sorted(listOfCloudEffectiveEmissivityPasses.keys()) :
+ emissSetNumber = emissName.split('_')[-1] # get the number off the end of the name
+ # if the numbering is at all different between the variables, this will fail
+ backgroundData[zCounter] = [listOfCloudEffectiveEmissivityPasses[emissName],
+ listOfOptionalLat[defaultValues['optLatBase'] + emissSetNumber],
+ listOfOptionalLon[defaultValues['optLonBase'] + emissSetNumber],
+ matplotlib.cm.gray, tempLevels]
+ zCounter = zCounter + 1
+
+ # if we're doing the optical depth land and ocean, add those to the background list too
+ if doOptionalDepthLandAndOcean :
+ tempLevels = _make_range_with_masked_arrays(listOfOpticalDepthLandAndOcean.values(), LEVELS_FOR_BACKGROUND_PLOTS)
+ for aodName in sorted(listOfOpticalDepthLandAndOcean.keys()) :
+ aodNumber = aodName.split('_')[-1] # get the number off the end of the name
+ # if the numbering is at all different between the variables, this will fail
+ backgroundData[zCounter] = [listOfOpticalDepthLandAndOcean[aodName],
+ listOfOptionalLat[defaultValues['optLatBase'] + aodNumber],
+ listOfOptionalLon[defaultValues['optLonBase'] + aodNumber],
+ matplotlib.cm.jet, tempLevels]
+ zCounter = zCounter + 1
+
# make the plot
LOG.info("Creating trajectory plot for time " + str(float(currentTime)) + ".")
tempFigure = _create_imapp_figure (initAODFlat, initLonData, initLatData,
pressureData=thisFramePressures, pressLongitudeData=thisFramePressLon, pressLatitudeData=thisFramePressLat,
- baseMapInstance=basemapObject, latRange=latRange, lonRange=lonRange)
+ baseMapInstance=basemapObject, latRange=latRange, lonRange=lonRange,
+ windsDataU=currentWindsU, windsDataV=currentWindsV,
+ windsDataLon=optionalDataWindsLongitude, windsDataLat=optionalDataWindsLatitude,
+ figureTitle="MODIS AOD & AOD Trajectories on " + str(currentTime) + "Z",
+ backgroundDataSets=backgroundData)
+ # TODO, the figure title needs to have the dates in it and the day may roll over, so current time can't be put in directly
# save the plot to disk
LOG.info("Saving plot to disk.")