# Copyright (C) 2012 VT SuperDARN Lab
# Full license can be found in LICENSE.txt
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
.. module:: fan
:synopsis: A module generating fan plots
.. moduleauthor:: AJ, 20130218
***************************
**Module**: pydarn.plotting.fan
***************************
**Functions**:
* :func:`pydarn.plotting.fan.plotFan`
* :func:`pydarn.plotting.fan.overlayFan`
"""
import pydarn,numpy,math,matplotlib,calendar,datetime,utils,pylab
import matplotlib.pyplot as plot
import matplotlib.lines as lines
from matplotlib.ticker import MultipleLocator
import matplotlib.patches as patches
from matplotlib.collections import PolyCollection,LineCollection
from mpl_toolkits.basemap import Basemap, pyproj
from utils.timeUtils import *
from pydarn.sdio.radDataRead import *
from matplotlib.figure import Figure
import matplotlib.cm as cm
from matplotlib.backends.backend_agg import FigureCanvasAgg
[docs]def plotFan(sTime,rad,interval=60,fileType='fitex',param='velocity',filtered=False ,\
scale=[],channel='a',coords='geo',colors='lasse',gsct=False,fov=True,edgeColors='face',lowGray=False,fill=True,\
velscl=1000.,legend=True,overlayPoes=False,poesparam='ted',poesMin=-3.,poesMax=0.5, \
poesLabel=r"Total Log Energy Flux [ergs cm$^{-2}$ s$^{-1}$]",overlayBnd=False, \
show=True,png=False,pdf=False,dpi=500):
"""A function to make a fan plot
**Args**:
* **sTime** (`datetime <http://tinyurl.com/bl352yx>`_): the start time you want to plot
* **rad** (list): a list of 3 letter radar codes, e.g. ['bks'], e.g. ['bks','wal','gbr']
* **[interval]** (int): the the time period to be plotted, in seconds. default = 60
* **[fileType]** (str): the file type to plot, valid inputs are 'fitex','fitacf', 'lmfit'. default = 'fitex'
* **[param]** (str): the parameter to be plotted, valid inputs are 'velocity', 'power', 'width', 'elevation', 'phi0'. default = 'velocity'
* **[filtered]** (boolean): a flag indicating whether the data should be boxcar filtered. default = False
* **[scale]** (list): the min and max values of the color scale, i.e. [min,max]. If this is set to [], then default values will be used
* **[channel] (char)**: the channel for which to plot data. default = 'a'
* **[coords]** (str): the coordinate system to use, valid inputs are 'geo', 'mag'. default = 'geo'
* **[colors]** (str): the color map to use, valid inputs are 'lasse', 'aj'. default = 'lasse'
* **[gsct]** (boolean): a flag indicating whether to plot ground scatter as gray. default = False
* **[fov]** (boolean): a flag indicating whether to overplot the radar fields of view. default = True
* **[edgeColors]** (str): edge colors of the polygons, default = 'face'
* **[lowGray]** (boolean): a flag indicating whether to plot low velocities in gray. default = False
* **[fill]** (boolean): a flag indicating whether to plot filled or point RB cells. default = True
* **[velscl]** (float): the velocity to use as baseline for velocity vector length, only applicable if fill = 0. default = 1000
* **[legend]** (boolean): a flag indicating whether to plot the legend, only applicable if fill = 0. default = True
* **[overlayPoes]** (boolean): a flag indicating whether to overlay poes data. default = False
* **[poesparam]** (str): the poes parameter to plot. default = 'ted'. available params can be found in :class:`gme.sat.poes.poesRec`
* **[poesMin]** (float): the min value for the poes data color scale. default = -3.
* **[poesMax]** (float): the max value for the poes data color scale. default = 0.5
* **[poesLabel]** (str): the label for the poes color bar. default = r"Total Log Energy Flux [ergs cm$^{-2}$ s$^{-1}$]"
* **[overlayBnd]** (boolean): a flag indicating whether to plot an auroral boundary determined from fitting poes data. default = False
* **[show]** (boolean): a flag indicating whether to display the figure on the screen. This can cause problems over ssh. default = True
* **[pdf]** (boolean): a flag indicating whether to output to a pdf file. default = False. WARNING: saving as pdf is slow
* **[png]** (boolean): a flag indicating whether to output to a png file. default = False
* **[dpi]** (int): dots per inch if saving as png. default = 300
**Returns**:
* Nothing
**Example**:
::
import datetime as dt
pydarn.plotting.fan.plotFan(dt.datetime(2013,3,16,16,30),['fhe','fhw'],param='power',gsct=True)
Written by AJ 20121004
"""
import datetime as dt, gme, pickle
from matplotlib.backends.backend_pdf import PdfPages
import models.aacgm as aacgm, os, copy
tt = dt.datetime.now()
#check the inputs
assert(isinstance(sTime,dt.datetime)),'error, sTime must be a datetime object'
assert(isinstance(rad,list)),"error, rad must be a list, eg ['bks'] or ['bks','fhe']"
for r in rad:
assert(isinstance(r,str) and len(r) == 3),'error, elements of rad list must be 3 letter strings'
assert(coords == 'geo' or coords == 'mag'),"error, coords must be one of 'geo' or 'mag'"
assert(param == 'velocity' or param == 'power' or param == 'width' or \
param == 'elevation' or param == 'phi0'), \
"error, allowable params are 'velocity','power','width','elevation','phi0'"
assert(scale == [] or len(scale)==2), \
'error, if present, scales must have 2 elements'
assert(colors == 'lasse' or colors == 'aj'),"error, valid inputs for color are 'lasse' and 'aj'"
if(scale == []):
if(param == 'velocity'): scale=[-200,200]
elif(param == 'power'): scale=[0,30]
elif(param == 'width'): scale=[0,150]
elif(param == 'elevation'): scale=[0,50]
elif(param == 'phi0'): scale=[-numpy.pi,numpy.pi]
fbase = sTime.strftime("%Y%m%d")
cmap,norm,bounds = utils.plotUtils.genCmap(param,scale,colors=colors,lowGray=lowGray)
#open the data files
myFiles = []
for r in rad:
f = radDataOpen(sTime,r,sTime+datetime.timedelta(seconds=interval),fileType=fileType,filtered=filtered,channel=channel)
if(f != None): myFiles.append(f)
assert(myFiles != []),'error, no data available for this period'
xmin,ymin,xmax,ymax = 1e16,1e16,-1e16,-1e16
allBeams = [''] * len(myFiles)
sites,fovs,oldCpids,lonFull,latFull=[],[],[],[],[]
#go through all open files
for i in range(len(myFiles)):
#read until we reach start time
allBeams[i] = radDataReadRec(myFiles[i])
while(allBeams[i].time < sTime and allBeams[i] != None):
allBeams[i] = radDataReadRec(myFiles[i])
#check that the file has data in the target interval
if(allBeams[i] == None):
myFiles[i].close()
myFiles[i] = None
continue
#get to field of view coords in order to determine map limits
t=allBeams[i].time
site = pydarn.radar.site(radId=allBeams[i].stid,dt=t)
sites.append(site)
if(coords == 'geo'):
latFull.append(site.geolat)
lonFull.append(site.geolon)
elif(coords == 'mag'):
x = aacgm.aacgmConv(site.geolat,site.geolon,0.,0)
latFull.append(x[0])
lonFull.append(x[1])
myFov = pydarn.radar.radFov.fov(site=site,rsep=allBeams[i].prm.rsep,\
ngates=allBeams[i].prm.nrang+1,nbeams=site.maxbeam,coords=coords)
fovs.append(myFov)
for b in range(0,site.maxbeam+1):
for k in range(0,allBeams[i].prm.nrang+1):
lonFull.append(myFov.lonFull[b][k])
latFull.append(myFov.latFull[b][k])
oldCpids.append(allBeams[i].cp)
#do some stuff in map projection coords to get necessary width and height of map
t1=dt.datetime.now()
lonFull,latFull = (numpy.array(lonFull)+360.)%360.,numpy.array(latFull)
tmpmap = Basemap(projection='npstere', boundinglat=20,lat_0=90, lon_0=numpy.mean(lonFull))
x,y = tmpmap(lonFull,latFull)
minx = x.min()
miny = y.min()
maxx = x.max()
maxy = y.max()
width = (maxx-minx)
height = (maxy-miny)
cx = minx + width/2.
cy = miny + height/2.
lon_0,lat_0 = tmpmap(cx, cy, inverse=True)
dist = width/50.
cTime = sTime
# if show:
# myFig = plot.figure(figsize=(12,8))
# else:
# myFig = Figure(figsize=(12,8))
myFig = plot.figure(figsize=(12,8))
#draw the actual map we want
myMap = Basemap(projection='stere',width=width,height=height,lon_0=numpy.mean(lonFull),lat_0=lat_0)
myMap.drawparallels(numpy.arange(-80.,81.,10.),labels=[1,0,0,0])
myMap.drawmeridians(numpy.arange(-180.,181.,20.),labels=[0,0,0,1])
if(coords == 'geo'):
myMap.drawcoastlines(linewidth=0.5,color='k')
myMap.drawmapboundary(fill_color='w')
myMap.fillcontinents(color='w', lake_color='w')
#overlay fields of view, if desired
if(fov == 1):
for r in rad:
pydarn.plotting.overlayRadar(myMap, codes=r, dateTime=sTime)
pydarn.plotting.overlayFov(myMap, codes=r, dateTime=sTime)
print dt.datetime.now()-t1
#manually draw the legend
if((not fill) and legend):
#draw the box
y = [myMap.urcrnry*.82,myMap.urcrnry*.99]
x = [myMap.urcrnrx*.86,myMap.urcrnrx*.99]
verts = [x[0],y[0]],[x[0],y[1]],[x[1],y[1]],[x[1],y[0]]
poly = patches.Polygon(verts,fc='w',ec='k',zorder=11)
myFig.gca().add_patch(poly)
labs = ['5 dB','15 dB','25 dB','35 dB','gs','1000 m/s']
pts = [5,15,25,35]
#plot the icons and labels
for w in range(6):
myFig.gca().text(x[0]+.35*(x[1]-x[0]),y[1]*(.98-w*.025),labs[w],zorder=15,color='k',size=8,va='center')
xctr = x[0]+.175*(x[1]-x[0])
if(w < 4):
myFig.scatter(xctr,y[1]*(.98-w*.025),s=.1*pts[w],zorder=15,marker='o',linewidths=.5,\
edgecolor='face',facecolor='k')
elif(w == 4):
myFig.scatter(xctr,y[1]*(.98-w*.025),s=.1*35.,zorder=15,marker='o',\
linewidths=.5,edgecolor='k',facecolor='w')
elif(w == 5):
y=LineCollection(numpy.array([((xctr-dist/2.,y[1]*(.98-w*.025)),(xctr+dist/2.,y[1]*(.98-w*.025)))]),linewidths=.5,zorder=15,color='k')
myFig.gca().add_collection(y)
# pickle.dump(myFig,open('map.pickle','wb'),-1)
bbox = myFig.gca().get_axes().get_position()
#now, loop through desired time interval
tz = dt.datetime.now()
cols = []
bndTime = sTime + datetime.timedelta(seconds=interval)
ft = 'None'
#go though all files
for i in range(len(myFiles)):
scans = [[] for j in range(len(myFiles))]
#check that we have good data at this time
if(myFiles[i] == None or allBeams[i] == None): continue
ft = allBeams[i].fType
#until we reach the end of the time window
while(allBeams[i] != None and allBeams[i].time < bndTime):
scans[i].append(allBeams[i])
#read the next record
allBeams[i] = radDataReadRec(myFiles[i])
intensities, pcoll = overlayFan(scans[i],myMap,myFig,param,coords,gsct=gsct,site=sites[i],fov=fovs[i], fill=fill,velscl=velscl,dist=dist,cmap=cmap,norm=norm)
cbar = myFig.colorbar(pcoll,orientation='vertical',shrink=.65,fraction=.1,drawedges=True)
l = []
#define the colorbar labels
for i in range(0,len(bounds)):
if(param == 'phi0'):
ln = 4
if(bounds[i] == 0): ln = 3
elif(bounds[i] < 0): ln = 5
l.append(str(bounds[i])[:ln])
continue
if((i == 0 and param == 'velocity') or i == len(bounds)-1):
l.append(' ')
continue
l.append(str(int(bounds[i])))
cbar.ax.set_yticklabels(l)
cbar.ax.tick_params(axis='y',direction='out')
#set colorbar ticklabel size
for ti in cbar.ax.get_yticklabels():
ti.set_fontsize(12)
if(param == 'velocity'):
cbar.set_label('Velocity [m/s]',size=14)
cbar.extend='max'
if(param == 'grid'): cbar.set_label('Velocity [m/s]',size=14)
if(param == 'power'): cbar.set_label('Power [dB]',size=14)
if(param == 'width'): cbar.set_label('Spec Wid [m/s]',size=14)
if(param == 'elevation'): cbar.set_label('Elev [deg]',size=14)
if(param == 'phi0'): cbar.set_label('Phi0 [rad]',size=14)
#myFig.gca().set_rasterized(True)
#label the plot
tx1 = myFig.text((bbox.x0+bbox.x1)/2.,bbox.y1+.02,cTime.strftime('%Y/%m/%d'),ha='center',size=14,weight=550)
tx2 = myFig.text(bbox.x1,bbox.y1+.02,cTime.strftime('%H:%M - ')+\
bndTime.strftime('%H:%M '),ha='right',size=13,weight=550)
tx3 = myFig.text(bbox.x0,bbox.y1+.02,'['+ft+']',ha='left',size=13,weight=550)
if(overlayPoes):
pcols = gme.sat.poes.overlayPoesTed(myMap, myFig.gca(), cTime, param=poesparam, scMin=poesMin, scMax=poesMax)
if(pcols != None):
cols.append(pcols)
pTicks = numpy.linspace(poesMin,poesMax,8)
cbar = myFig.colorbar(pcols,ticks=pTicks,orientation='vertical',shrink=0.65,fraction=.1)
cbar.ax.set_yticklabels(pTicks)
cbar.set_label(poesLabel,size=14)
cbar.ax.tick_params(axis='y',direction='out')
#set colorbar ticklabel size
for ti in cbar.ax.get_yticklabels():
ti.set_fontsize(12)
if(overlayBnd):
gme.sat.poes.overlayPoesBnd(myMap, myFig.gca(), cTime)
#handle the outputs
if png == True:
# if not show:
# canvas = FigureCanvasAgg(myFig)
myFig.savefig(sTime.strftime("%Y%m%d.%H%M.")+str(interval)+'.fan.png',dpi=dpi)
if pdf:
# if not show:
# canvas = FigureCanvasAgg(myFig)
myFig.savefig(sTime.strftime("%Y%m%d.%H%M.")+str(interval)+'.fan.pdf')
if show:
myFig.show()
[docs]def overlayFan(myData,myMap,myFig,param,coords='geo',gsct=0,site=None,\
fov=None,gs_flg=[],fill=True,velscl=1000.,dist=1000.,
cmap=None,norm=None,alpha=1):
"""A function of overlay radar scan data on a map
**Args**:
* **myData (:class:`pydarn.sdio.radDataTypes.scanData` or :class:`pydarn.sdio.radDataTypes.beamData` or list of :class:`pydarn.sdio.radDataTypes.beamData` objects)**: a radar beam object, a radar scanData object, or simply a list of radar beams
* **myMap**: the map we are plotting on
* **[param]**: the parameter we are plotting
* **[coords]**: the coordinates we are plotting in
* **[param]**: the parameter to be plotted, valid inputs are 'velocity', 'power', 'width', 'elevation', 'phi0'. default = 'velocity
* **[gsct]**: a flag indicating whether we are distinguishing ground scatter. default = 0
* **[intensities]**: a list of intensities (used for colorbar)
* **[fov]**: a radar fov object
* **[gs_flg]**: a list of gs flags, 1 per range gate
* **[fill]**: a flag indicating whether to plot filled or point RB cells. default = True
* **[velscl]**: the velocity to use as baseline for velocity vector length, only applicable if fill = 0. default = 1000
* **[lines]**: an array to have the endpoints of velocity vectors. only applicable if fill = 0. default = []
* **[dist]**: the length in map projection coords of a velscl length velocity vector. default = 1000. km
**OUTPUTS**:
NONE
**EXAMPLE**:
::
overlayFan(aBeam,myMap,param,coords,gsct=gsct,site=sites[i],fov=fovs[i],\
verts=verts,intensities=intensities,gs_flg=gs_flg)
Written by AJ 20121004
"""
if(site == None):
site = pydarn.radar.site(radId=myData[0].stid, dt=myData[0].time)
if(fov == None):
fov = pydarn.radar.radFov.fov(site=site,rsep=myData[0].prm.rsep,\
ngates=myData[0].prm.nrang+1,nbeams= site.maxbeam,coords=coords)
if(isinstance(myData,pydarn.sdio.beamData)): myData = [myData]
gs_flg,lines = [],[]
if fill: verts,intensities = [],[]
else: verts,intensities = [[],[]],[[],[]]
#loop through gates with scatter
for myBeam in myData:
for k in range(0,len(myBeam.fit.slist)):
if myBeam.fit.slist[k] not in fov.gates: continue
r = myBeam.fit.slist[k]
if fill:
x1,y1 = myMap(fov.lonFull[myBeam.bmnum,r],fov.latFull[myBeam.bmnum,r])
x2,y2 = myMap(fov.lonFull[myBeam.bmnum,r+1],fov.latFull[myBeam.bmnum,r+1])
x3,y3 = myMap(fov.lonFull[myBeam.bmnum+1,r+1],fov.latFull[myBeam.bmnum+1,r+1])
x4,y4 = myMap(fov.lonFull[myBeam.bmnum+1,r],fov.latFull[myBeam.bmnum+1,r])
#save the polygon vertices
verts.append(((x1,y1),(x2,y2),(x3,y3),(x4,y4),(x1,y1)))
#save the param to use as a color scale
if(param == 'velocity'): intensities.append(myBeam.fit.v[k])
elif(param == 'power'): intensities.append(myBeam.fit.p_l[k])
elif(param == 'width'): intensities.append(myBeam.fit.w_l[k])
elif(param == 'elevation' and myBeam.prm.xcf): intensities.append(myBeam.fit.elv[k])
elif(param == 'phi0' and myBeam.prm.xcf): intensities.append(myBeam.fit.phi0[k])
else:
x1,y1 = myMap(fov.lonCenter[myBeam.bmnum,r],fov.latCenter[myBeam.bmnum,r])
verts[0].append(x1)
verts[1].append(y1)
x2,y2 = myMap(fov.lonCenter[myBeam.bmnum,r+1],fov.latCenter[myBeam.bmnum,r+1])
theta = math.atan2(y2-y1,x2-x1)
x2,y2 = x1+myBeam.fit.v[k]/velscl*(-1.0)*math.cos(theta)*dist,y1+myBeam.fit.v[k]/velscl*(-1.0)*math.sin(theta)*dist
lines.append(((x1,y1),(x2,y2)))
#save the param to use as a color scale
if(param == 'velocity'): intensities[0].append(myBeam.fit.v[k])
elif(param == 'power'): intensities[0].append(myBeam.fit.p_l[k])
elif(param == 'width'): intensities[0].append(myBeam.fit.w_l[k])
elif(param == 'elevation' and myBeam.prm.xcf): intensities[0].append(myBeam.fit.elv[k])
elif(param == 'phi0' and myBeam.prm.xcf): intensities[0].append(myBeam.fit.phi0[k])
if(myBeam.fit.p_l[k] > 0): intensities[1].append(myBeam.fit.p_l[k])
else: intensities[1].append(0.)
if(gsct): gs_flg.append(myBeam.fit.gflg[k])
#do the actual overlay
if(fill):
#if we have data
if(verts != []):
if(gsct == 0):
inx = numpy.arange(len(verts))
else:
inx = numpy.where(numpy.array(gs_flg)==0)
x = PolyCollection(numpy.array(verts)[numpy.where(numpy.array(gs_flg)==1)],
facecolors='.3',linewidths=0,zorder=5,alpha=alpha)
myFig.gca().add_collection(x, autolim=True)
pcoll = PolyCollection(numpy.array(verts)[inx],
edgecolors='face',linewidths=0,closed=False,zorder=4,
alpha=alpha,cmap=cmap,norm=norm)
#set color array to intensities
pcoll.set_array(numpy.array(intensities)[inx])
myFig.gca().add_collection(pcoll, autolim=True)
return intensities,pcoll
else:
#if we have data
if(verts != [[],[]]):
if(gsct == 0):
inx = numpy.arange(len(verts[0]))
else:
inx = numpy.where(numpy.array(gs_flg)==0)
#plot the ground scatter as open circles
x = myFig.scatter(numpy.array(verts[0])[numpy.where(numpy.array(gs_flg)==1)],\
numpy.array(verts[1])[numpy.where(numpy.array(gs_flg)==1)],\
s=.1*numpy.array(intensities[1])[numpy.where(numpy.array(gs_flg)==1)],\
zorder=5,marker='o',linewidths=.5,facecolors='w',edgecolors='k')
myFig.gca().add_collection(x, autolim=True)
#plot the i-s as filled circles
ccoll = myFig.gca().scatter(numpy.array(verts[0])[inx],numpy.array(verts[1])[inx], \
s=.1*numpy.array(intensities[1])[inx],zorder=10,marker='o',linewidths=.5, \
edgecolors='face',cmap=cmap,norm=norm)
#set color array to intensities
ccoll.set_array(numpy.array(intensities[0])[inx])
myFig.gca().add_collection(ccoll)
#plot the velocity vectors
lcoll = LineCollection(numpy.array(lines)[inx],linewidths=.5,zorder=12,cmap=cmap,norm=norm)
lcoll.set_array(numpy.array(intensities[0])[inx])
myFig.gca().add_collection(lcoll)
return intensities,lcoll