[cig-commits] r11779 - cs/portal/trunk/seismo/SeismoWebPortal
leif at geodynamics.org
leif at geodynamics.org
Tue Apr 8 22:16:07 PDT 2008
Author: leif
Date: 2008-04-08 22:16:07 -0700 (Tue, 08 Apr 2008)
New Revision: 11779
Modified:
cs/portal/trunk/seismo/SeismoWebPortal/gmt.py
cs/portal/trunk/seismo/SeismoWebPortal/views.py
Log:
Incorporated Vala's plot_wig_map.
Modified: cs/portal/trunk/seismo/SeismoWebPortal/gmt.py
===================================================================
--- cs/portal/trunk/seismo/SeismoWebPortal/gmt.py 2008-04-09 01:06:43 UTC (rev 11778)
+++ cs/portal/trunk/seismo/SeismoWebPortal/gmt.py 2008-04-09 05:16:07 UTC (rev 11779)
@@ -164,6 +164,432 @@
return jpg
+def wiggleMap(cmt,
+ stationList,
+ opt_a=None,
+ opt_A=None,
+ opt_B=None,
+ opt_c=None,
+ opt_D=None,
+ opt_G=None,
+ opt_l=None,
+ opt_L=None,
+ opt_e=None, opt_m=None,
+ opt_M=None,
+ opt_R=None,
+ opt_s=None,
+ opt_S=None,
+ opt_t=None,
+ opt_W=None,
+ opt_V=None,
+ ):
+ # Based on:
+ # plot_wig_map.pl
+ # Vala Hjorleifsdottir
+ # Feb 27 2003
+
+ # align_phase -- one of o,p,s,r,l,number, aligns phase with station loc.
+ align = None
+ if opt_a is not None:
+ if opt_a in "opsrl":
+ align = opt_a
+ else:
+ try:
+ align = float(opt_a)
+ except ValueError:
+ assert False, "Align on o,p,s,r,l or number, not: %r" % opt_a
+
+ # plot antipode [useful for global map]
+ if opt_A and not opt_R:
+ plot_antipode = True
+ elif not opt_A and not opt_R:
+ plot_antipode = False
+ arc_min = -1
+ arc_max = 90
+
+ # plot both sides of globe [also for global maps]
+ if opt_B:
+ nr_of_plots = 2
+ else:
+ nr_of_plots = 1
+
+ # cut0/cut1 -- plot data from cut1 to cut2 relative to align_phase [no cut]
+ if opt_c:
+ cut_option = "-C%s/%s" % opt_c
+ else:
+ cut_option = ""
+
+ # xshift/yshift -- shifts trace by these amounts of MEASURE_UNIT [no shift]
+ if opt_D:
+ shift_option = "-D%s" % opt_D
+ else:
+ shift_option = "-D0/0"
+
+ # plot great circle path
+ if opt_G:
+ plot_greatcircle = True
+ else:
+ plot_greatcircle = False
+
+ # tscale -- seconds/inch, controls length of trace [3600]
+ if opt_l:
+ tscale = opt_l
+ else:
+ tscale = "3600"
+
+ # Landscape
+ if opt_L:
+ landscape = True
+ else:
+ landscape = False
+
+ # Lat/Lon -- give either -e or -m to determine center of map
+ if not opt_m and not opt_e:
+ assert False, "You have to specify either -e Lat/Lon or -m CMTFILE"
+ if opt_e:
+ plot_mecha = False
+ # CMT-file -- CMT-file in HVD format, plots mechanism on map
+ if opt_m:
+ plot_mecha = True
+
+ # scale -- scale for seismograms, type pssac2 for options
+ if opt_M:
+ wig_scale = opt_M
+ else:
+ wig_scale = "N/A"
+
+ # lon1/lon2/lat1/lat2/tick_x/tick_y/width/x/y - regional map. Width of map and x and y of lower left corner [in inches].
+ if opt_R:
+ regional = True
+ min_lon, max_lon, min_lat, max_lat, tick_x, tick_y, scale, x, y = opt_R
+ plot_antipode = False
+ arc_min = 0
+ arc_max = 180
+ else:
+ regional = False
+ min_lon, max_lon, min_lat, max_lat, tick_x, tick_y, scale, x, y = None, None, None, None, None, None, None, None, None
+
+ # plot stations
+ if opt_s:
+ plot_stations = True
+ else:
+ plot_stations = False
+
+ # syndir,synsuff -- plot synthetics SYNDIR/*synsuff correpsonding to datafiles
+ if opt_S:
+ syn_dir, syn_suff = opt_S
+ assert os.path.isdir(syn_dir), "No such directory: %s" % syn_dir
+ plot_syn = True
+ else:
+ syn_dir, syn_suff = None, None
+ plot_syn = False
+
+ # title -- title of plot [no title]
+ if opt_t:
+ title = opt_t
+ else:
+ title = ""
+
+ # plot wiggles (give datafiles at the end)
+ if opt_W:
+ plot_wiggles = True
+ else:
+ plot_wiggles = False
+
+ assert not landscape or regional, "You can only use landscape for regional plots"
+ assert opt_M or not plot_wiggles, "Specify scale for wiggles"
+ if plot_wiggles and plot_syn and not '/' in opt_M:
+ assert False, "Don't use relative scale when plotting both data and synthetics -- change -M option"
+ assert plot_wiggles or not plot_syn, "Have to specify datafiles to pick corresponding synthetics"
+
+ # verbose
+ if opt_V:
+ if not opt_l and plot_wiggles:
+ print "Using default timescale, 1hour/inch"
+ print "Aligning on: %s" % align
+ print "Plotting antipode? %s" % plot_antipode
+ print "Cutting from/to sec: %s" % opt_c
+ #print "Moment tensor file: %s" % cmt_file
+ print "Scale for wiggles: %s" % wig_scale
+ print "Name of outputfile: %s" % PS_OUT
+ print "Plotting stations? %s" % plot_stations
+ print "Plotting synthetics? %s" % plot_syn
+ print "Title of plot: %s" % title
+ print "Plotting wiggles: %s" % plot_wiggles
+ print "Scale for wiggles: %s" % wig_scale
+ V_option=" -V "
+ else:
+ V_option = ""
+
+
+ if plot_mecha:
+ event_lat = cmt.latitude
+ event_lon = cmt.longitude
+ else:
+ event_lat, event_lon = opt_e
+
+ args = [
+ align, plot_antipode,
+ arc_min, arc_max,
+ cut_option, shift_option,
+ plot_greatcircle, tscale, landscape, plot_mecha,
+ wig_scale, regional,
+ min_lon, max_lon, min_lat, max_lat, tick_x, tick_y, scale, x, y,
+ plot_stations,
+ plot_syn, syn_dir, syn_suff,
+ title,
+ plot_wiggles,
+ V_option,
+ event_lat, event_lon,
+ opt_R,
+ ]
+
+ owd = os.getcwd()
+ workdir = mkdtemp()
+ os.chdir(workdir)
+ os.putenv("LD_LIBRARY_PATH", "/home/leif/opt/netCDF/lib")
+
+ PS_OUT = join(workdir, "map.ps")
+
+ plotWiggleMap(1, nr_of_plots, cmt, stationList, PS_OUT, *args)
+ if nr_of_plots == 2:
+ plotWiggleMap(2, nr_of_plots, cmt, stationList, PS_OUT, *args)
+
+ os.chdir(owd)
+
+ # Convert the PostScript to a JPG file.
+ jpg = join(workdir, "map.jpg")
+ spawn("/usr/bin/convert", '-trim', '+repage', PS_OUT, jpg)
+
+ jpg = open(jpg, 'r')
+ shutil.rmtree(workdir)
+
+ return jpg
+
+
+def plotWiggleMap(plot_nr, nr_of_plots,
+ cmt, stationList, PS_OUT,
+ # *args
+ align, plot_antipode,
+ arc_min, arc_max,
+ cut_option, shift_option,
+ plot_greatcircle, tscale, landscape, plot_mecha,
+ wig_scale, regional,
+ min_lon, max_lon, min_lat, max_lat, tick_x, tick_y, scale, x, y,
+ plot_stations,
+ plot_syn, syn_dir, syn_suff,
+ title,
+ plot_wiggles,
+ V_option,
+ event_lat, event_lon,
+ opt_R,
+ ):
+
+ #########################
+ #Deciding on plot options
+ #########################
+
+ #@data_files=@data_files_original;
+
+ if nr_of_plots > 1:
+ landscape = True
+ if plot_nr == 1:
+ plot_antipode = False
+ elif plot_nr == 2:
+ plot_antipode = True
+ else:
+ assert False, "confusing plot_nr: $plot_nr - something broke in script"
+
+ if plot_antipode:
+ arc_min = 90
+ arc_max = 180
+ event_lat = -event_lat
+ if event_lon > 0:
+ event_lon = event_lon-180
+ else:
+ event_lon = event_lon+180
+ elif not regional:
+ arc_min = -1
+ arc_max = 90
+
+ if regional:
+ x=x+1
+ y=y+2
+ bounds="-R%s/%s/%s/%s" % (min_lon, max_lon, min_lat, max_lat)
+ proj = "-JM%s" % scale
+ tick = "-B%s/%sSWen" % (tick_x, tick_y)
+ if max_lon-min_lon > 30:
+ resolution="-Dl -A1000"
+ elif max_lon-min_lon<=30 and max_lon-min_lon>5:
+ resolution="-Di -A100"
+ elif max_lon-min_lon<=5 and max_lon-min_lon>1:
+ resolution="-Dh -A10"
+ else:
+ resolution="-Df -A1"
+ else:
+ if nr_of_plots==1:
+ scale = 7.0
+ x=0.7
+ y=1.8
+ else:
+ scale=4.9
+ if plot_nr==1:
+ x=0.5
+ y=1.8
+ else:
+ x=5.6
+ y=1.8
+ proj = "-JA%s/%s/%s" % (event_lon, event_lat, scale)
+ bounds = "-R-180/180/-90/90"
+ tick = "-B60/60"
+ resolution="-Dl -A1000"
+
+ if landscape:
+ p_option=""
+ # l="-l4.8/-0.2/1/0.1/8"
+ if nr_of_plots==1:
+ lx=5.2-x
+ ly=1.2-y
+ time_loc=" 5.2 0.8 "
+ tit_loc=" 5.2 7.2 "
+ else:
+ lx=5.5-x
+ ly=1.9-y
+ time_loc=" 5.5 1.5 "
+ tit_loc=" 5.5 6.7 "
+ l="-l%s/%s/1/0.1/8" % (lx, ly)
+ else:
+ p_option="-P"
+ lx=4.2-x
+ ly=1.6-y
+ l="-l%s/%s/1/0.1/8" % (lx, ly)
+ time_loc=" 4.2 1.2 "
+ tit_loc=" 4.2 9.2 "
+
+ scale2 = 1.10*scale/3
+ scale3 = 2.12*scale/3
+
+ bin = gmt_bin
+
+ # set GMT default commands
+ os.system("%(bin)s/gmtset MEASURE_UNIT inch DEGREE_FORMAT 3 BASEMAP_TYPE plain DOTS_PR_INCH 600 PAPER_MEDIA letter+" % locals())
+
+ if isinstance(align, float):
+ e_label = "-En%f" % align
+ time_label = "Time-distance/%f (km/s) [s]" % align
+ else:
+ alignMap = {
+ 'o': ("-Ent-3", "Time (s) from event"),
+ 'p': ("-Ent1", "Time (s) (aligned on p)"),
+ 's': ("-Ent2", "Time (s) (aligned on s)"),
+ 'r': ("-En3.8", "Time-distance/3.8(km/s) (s)"),
+ 'l': ("-En4.4", "Time-distance/4.4(km/s) (s)"),
+ None: (None, None)
+ }
+ e_label, time_label = alignMap[align]
+
+ if plot_stations or plot_wiggles:
+ if not opt_R:
+ #@data_files=sac_files_prune("gcarc",$arc_min,$arc_max, at data_files);
+ pass
+ else:
+ #@data_files=sac_files_prune("stla",$min_lat,$max_lat, at data_files);
+ #@data_files=sac_files_prune("stlo",$min_lon,$max_lon, at data_files);
+ pass
+
+ #foreach $file (@data_files){chomp($file);}
+ #@stlat=sac_files_get_values('stla', at data_files);
+ #@stlon=sac_files_get_values('stlo', at data_files);
+ #@stname = sac_files_get_values('kstnm', at data_files);
+ #@stnetwk = sac_files_get_values('knetwk', at data_files);
+ #@stcomp=sac_files_get_values('kcmpnm', at data_files);
+ if plot_wiggles and plot_nr==1:
+ title="$title $stcomp[0]"
+
+ ### PLOTTING ###
+
+ ##
+ #Plotting basemap (with equal distance rings)-S150/200/255 -G200/255/150 -G240/240/240 -S255/255/255
+ ##
+ if plot_nr==1:
+ O_option=""
+ output=">"
+ else:
+ O_option="-O"
+ output=">>"
+
+ os.system("%(bin)s/pscoast %(bounds)s %(proj)s %(tick)s %(resolution)s -W1/0/0/255 -S150/200/255 -G200/255/150 -K %(O_option)s %(V_option)s -X%(x)s -Y%(y)s %(p_option)s %(output)s %(PS_OUT)s" % locals())
+ if not regional:
+ inFile("%s %s" % (event_lon, event_lat))
+ os.system("%(bin)s/psxy %(bounds)s %(proj)s -Sc%(scale2)s -W2/0 -K -O %(V_option)s >>%(PS_OUT)s <inFile" % locals())
+
+ ##
+ #Plotting the event_location:
+ ##
+
+ if not plot_antipode and plot_mecha:
+ inFile("%f %f 0 %f %f %f %f %f %f 1 0 0" % (
+ event_lon, event_lat, cmt.Mrr, cmt.Mtt, cmt.Mpp, cmt.Mrt, cmt.Mrp, cmt.Mtp))
+ os.system("%(bin)s/psmeca %(bounds)s %(proj)s -Sm0.3 -L2 -G255/0/0 -T0 -K -O %(V_option)s >> %(PS_OUT)s <inFile" % locals())
+ else:
+ inFile("%s %s" % (event_lon, event_lat))
+ os.system("%(bin)s/psxy %(bounds)s %(proj)s -Sc0.1 -G255/255/0 -W2/0 -K -O %(V_option)s >>%(PS_OUT)s <inFile" % locals())
+
+ if plot_stations:
+ #Plotting the stations:
+ #open(PSXYst,"|psxy $bounds $proj -St0.1 -G255/0/0 -W2/0 -K -O $V_option >> %(PS_OUT)s");
+ #for($i = 0; $i < @data_files; $i++){print PSXYst "$stlon[$i] $stlat[$i] 8 0 4 CB \n";}
+ #close(PSXYst);
+ #Plotting the station-names:
+ #open(PSXYstn,"| pstext $bounds $proj -D0/0.05 -N -K -O -V_option >> %(PS_OUT)s");
+ #for($i = 0; $i < @data_files; $i++){print PSXYstn "$stlon[$i] $stlat[$i] 8 0 4 CB $stname[$i] \n";}
+ #close(PSXYstn);
+ pass
+
+ if plot_greatcircle:
+ #open(PSXYgcarc,"|psxy $bounds $proj -W5/0 -K -M -O $V_option >> %(PS_OUT)s");
+ #for($i = 0; $i < @data_files; $i++){print PSXYgcarc ">\n$event_lon $event_lat\n$stlon[$i] $stlat[$i]\n";}
+ #close(PSXYgcarc);
+ pass
+
+ if plot_wiggles:
+ #print "Plotting data\n";
+ #print "$pssac2 $bounds $proj $shift_option $l -M$wig_scale -N $cut_option -W2p/0/0/0 -L$tscale $e_label -P -O -K -V \n";
+ #open(PSSAC2, "| $pssac2 $bounds $proj $shift_option $l -M$wig_scale -N $cut_option -W2p/0/0/0 -L$tscale $e_label -P -O -K -V >> %(PS_OUT)s;") or die "pssac2 failed";
+ #for($i = 0;$i<@data_files;$i++){print PSSAC2 "$data_files[$i] $stlon[$i] $stlat[$i] 0.5p/0/0/0\n";}
+ #close(PSSAC2);
+ pass
+
+ if plot_syn:
+ #print "Plotting synthetics \n";
+ #my(@syn_files);
+ #for($i=0; $i<@data_files; $i++){
+ # $syn_file="${syn_dir}/$stname[$i].$stnetwk[$i].$stcomp[$i]$syn_suff";
+ # if(-e $syn_file){
+ # push @syn_files, $syn_file;
+ # }
+ #}
+ #@sylat=sac_files_get_values('stla', at syn_files);
+ #@sylon=sac_files_get_values('stlo', at syn_files);
+ #open(PSSAC2, "| $pssac2 $bounds $proj $shift_option -M$wig_scale -N $cut_option -W2p/0/0/0 -L$tscale $e_label -P -O -K -V >> $PS_OUT") or die "pssac2 failed";
+ #for($i = 0;$i<@syn_files;$i++){print PSSAC2 "$syn_files[$i] $sylon[$i] $sylat[$i] 0.5p/255/0/0\n";}
+ #close(PSSAC2);
+ pass
+
+ if plot_wiggles or plot_syn:
+ #system("echo $time_loc 16 0 5 CM \"$time_label\" | pstext -R0/10/0/10 -Jx1 -X-$x -Y-$y -K -O $V_option>>%(PS_OUT)s");
+ #x=0
+ #y=0
+ pass
+
+
+ #Plotting title:
+ inFile('%s 24 0 5 CM "%s"' % (tit_loc, title))
+ os.system("%(bin)s/pstext -R0/10/0/10 -Jx1 -O %(V_option)s -X-%(x)s -Y-%(y)s >>%(PS_OUT)s <inFile" % locals())
+
+ return
+
+
def spawn(*argv):
status = os.spawnvp(os.P_WAIT, argv[0], argv)
if status != 0:
Modified: cs/portal/trunk/seismo/SeismoWebPortal/views.py
===================================================================
--- cs/portal/trunk/seismo/SeismoWebPortal/views.py 2008-04-09 01:06:43 UTC (rev 11778)
+++ cs/portal/trunk/seismo/SeismoWebPortal/views.py 2008-04-09 05:16:07 UTC (rev 11779)
@@ -1153,8 +1153,11 @@
return HttpResponseRedirect(postTrashRedirect)
-def img(src):
- return '<img src="%s">' % src
+def img(src, width=None, height=None):
+ wh = ''
+ if width is not None:
+ wh = 'width=%d height=%d' % (width, height)
+ return '<img src="%s"%s>' % (src, wh)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Events
@@ -1198,6 +1201,7 @@
navtree = gui.Directory("events", "Events", [
gui.Directory("seismograms", "Seismograms", [], url = url),
sources,
+ gui.File("map", "Map", url = url + "map/"),
gui.File("settings", "Settings", url = url + "settings/"),
gui.File("properties", "Properties", url = url + "properties/"),
])
@@ -1213,6 +1217,24 @@
if name == "sources":
return configSources(workspace, url, request, path, appWindow, desktop)
+ if name == "map":
+ if path:
+ name = path.pop(0)
+ if path or name != "map.jpg": raise Http404
+ from gmt import wiggleMap
+ import shutil
+ stream = wiggleMap(event.singleSource, workspace.stations, opt_m=True)
+ response = HttpResponse(mimetype='image/jpeg')
+ shutil.copyfileobj(stream, response)
+ stream.close()
+ return response
+
+ appWindow.path.append(navtree.index['map'])
+ map = gui.ChildWindow(url + "map/", "Map")
+ map.content = gui.StaticContent(img(url + "map/map.jpg"))
+ desktop.activeWindow.selectWindow(map)
+ return desktop
+
if name == "settings":
appWindow.path.append(navtree.index['settings'])
settings = gui.ChildWindow(url + "settings/", "Settings")
@@ -1800,7 +1822,7 @@
)
if name == "map":
- map.content = gui.StaticContent(img(url + "map.jpg"))
+ map.content = gui.StaticContent(img(url + "map.jpg", width=547, height=657))
desktop.activeWindow.selectWindow(map)
return desktop
More information about the cig-commits
mailing list