- Timestamp:
- Apr 21, 2010, 1:48:09 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
r7685 r7690 317 317 318 318 return bearing 319 320 321 ##322 # @brief Generate figures from quantities and gauges for each sww file.323 # @param plot_quantity ??324 # @param file_loc ??325 # @param report ??326 # @param reportname ??327 # @param surface ??328 # @param leg_label ??329 # @param f_list ??330 # @param gauges ??331 # @param locations ??332 # @param elev ??333 # @param gauge_index ??334 # @param production_dirs ??335 # @param time_min ??336 # @param time_max ??337 # @param time_unit ??338 # @param title_on ??339 # @param label_id ??340 # @param generate_fig ??341 # @param verbose??342 # @return (texfile2, elev_output)343 def generate_figures(plot_quantity, file_loc, report, reportname, surface,344 leg_label, f_list, gauges, locations, elev, gauge_index,345 production_dirs, time_min, time_max, time_unit,346 title_on, label_id, generate_fig, verbose):347 """ Generate figures based on required quantities and gauges for348 each sww file349 """350 from os import sep, altsep, getcwd, mkdir, access, F_OK, environ351 352 if generate_fig is True:353 from pylab import ion, hold, plot, axis, figure, legend, savefig, \354 xlabel, ylabel, title, close, subplot355 356 if surface is True:357 import pylab as p1358 import mpl3d.mplot3d as p3359 360 if report == True:361 texdir = getcwd()+sep+'report'+sep362 if access(texdir,F_OK) == 0:363 mkdir (texdir)364 if len(label_id) == 1:365 label_id1 = label_id[0].replace(sep,'')366 label_id2 = label_id1.replace('_','')367 texfile = texdir + reportname + '%s' % label_id2368 texfile2 = reportname + '%s' % label_id2369 texfilename = texfile + '.tex'370 fid = open(texfilename, 'w')371 372 if verbose: log.critical('Latex output printed to %s' % texfilename)373 else:374 texfile = texdir+reportname375 texfile2 = reportname376 texfilename = texfile + '.tex'377 fid = open(texfilename, 'w')378 379 if verbose: log.critical('Latex output printed to %s' % texfilename)380 else:381 texfile = ''382 texfile2 = ''383 384 p = len(f_list)385 n = []386 n0 = 0387 for i in range(len(f_list)):388 n.append(len(f_list[i].get_time()))389 if n[i] > n0: n0 = n[i]390 n0 = int(n0)391 m = len(locations)392 model_time = num.zeros((n0, m, p), num.float)393 stages = num.zeros((n0, m, p), num.float)394 elevations = num.zeros((n0, m, p), num.float)395 momenta = num.zeros((n0, m, p), num.float)396 xmom = num.zeros((n0, m, p), num.float)397 ymom = num.zeros((n0, m, p), num.float)398 speed = num.zeros((n0, m, p), num.float)399 bearings = num.zeros((n0, m, p), num.float)400 due_east = 90.0*num.ones((n0, 1), num.float)401 due_west = 270.0*num.ones((n0, 1), num.float)402 depths = num.zeros((n0, m, p), num.float)403 eastings = num.zeros((n0, m, p), num.float)404 min_stages = []405 max_stages = []406 min_momentums = []407 max_momentums = []408 max_xmomentums = []409 max_ymomentums = []410 min_xmomentums = []411 min_ymomentums = []412 max_speeds = []413 min_speeds = []414 max_depths = []415 model_time_plot3d = num.zeros((n0, m), num.float)416 stages_plot3d = num.zeros((n0, m), num.float)417 eastings_plot3d = num.zeros((n0, m),num.float)418 if time_unit is 'mins': scale = 60.0419 if time_unit is 'hours': scale = 3600.0420 421 ##### loop over each swwfile #####422 for j, f in enumerate(f_list):423 if verbose: log.critical('swwfile %d of %d' % (j, len(f_list)))424 425 starttime = f.starttime426 comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'427 fid_compare = open(comparefile, 'w')428 file0 = file_loc[j] + 'gauges_t0.csv'429 fid_0 = open(file0, 'w')430 431 ##### loop over each gauge #####432 for k in gauge_index:433 if verbose: log.critical('Gauge %d of %d' % (k, len(gauges)))434 435 g = gauges[k]436 min_stage = 10437 max_stage = 0438 max_momentum = max_xmomentum = max_ymomentum = 0439 min_momentum = min_xmomentum = min_ymomentum = 100440 max_speed = 0441 min_speed = 0442 max_depth = 0443 gaugeloc = str(locations[k])444 thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \445 + gaugeloc + '.csv'446 if j == 0:447 fid_out = open(thisfile, 'w')448 s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'449 fid_out.write(s)450 451 #### generate quantities #######452 for i, t in enumerate(f.get_time()):453 if time_min <= t <= time_max:454 w = f(t, point_id = k)[0]455 z = f(t, point_id = k)[1]456 uh = f(t, point_id = k)[2]457 vh = f(t, point_id = k)[3]458 depth = w-z459 m = sqrt(uh*uh + vh*vh)460 if depth < 0.001:461 vel = 0.0462 else:463 vel = m / (depth + 1.e-6/depth)464 bearing = calc_bearing(uh, vh)465 model_time[i,k,j] = (t + starttime)/scale #t/60.0466 stages[i,k,j] = w467 elevations[i,k,j] = z468 xmom[i,k,j] = uh469 ymom[i,k,j] = vh470 momenta[i,k,j] = m471 speed[i,k,j] = vel472 bearings[i,k,j] = bearing473 depths[i,k,j] = depth474 thisgauge = gauges[k]475 eastings[i,k,j] = thisgauge[0]476 s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \477 % (t, w, m, vel, z, uh, vh, bearing)478 fid_out.write(s)479 if t == 0:480 s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)481 fid_0.write(s)482 if t/60.0 <= 13920: tindex = i483 if w > max_stage: max_stage = w484 if w < min_stage: min_stage = w485 if m > max_momentum: max_momentum = m486 if m < min_momentum: min_momentum = m487 if uh > max_xmomentum: max_xmomentum = uh488 if vh > max_ymomentum: max_ymomentum = vh489 if uh < min_xmomentum: min_xmomentum = uh490 if vh < min_ymomentum: min_ymomentum = vh491 if vel > max_speed: max_speed = vel492 if vel < min_speed: min_speed = vel493 if z > 0 and depth > max_depth: max_depth = depth494 495 496 s = '%.2f, %.2f, %.2f, %.2f, %s\n' \497 % (max_stage, min_stage, z, thisgauge[0], leg_label[j])498 fid_compare.write(s)499 max_stages.append(max_stage)500 min_stages.append(min_stage)501 max_momentums.append(max_momentum)502 max_xmomentums.append(max_xmomentum)503 max_ymomentums.append(max_ymomentum)504 min_xmomentums.append(min_xmomentum)505 min_ymomentums.append(min_ymomentum)506 min_momentums.append(min_momentum)507 max_depths.append(max_depth)508 max_speeds.append(max_speed)509 min_speeds.append(min_speed)510 #### finished generating quantities for each swwfile #####511 512 model_time_plot3d[:,:] = model_time[:,:,j]513 stages_plot3d[:,:] = stages[:,:,j]514 eastings_plot3d[:,] = eastings[:,:,j]515 516 if surface is True:517 log.critical('Printing surface figure')518 for i in range(2):519 fig = p1.figure(10)520 ax = p3.Axes3D(fig)521 if len(gauges) > 80:522 ax.plot_surface(model_time[:,:,j],523 eastings[:,:,j],524 stages[:,:,j])525 else:526 ax.plot3D(num.ravel(eastings[:,:,j]),527 num.ravel(model_time[:,:,j]),528 num.ravel(stages[:,:,j]))529 ax.set_xlabel('time')530 ax.set_ylabel('x')531 ax.set_zlabel('stage')532 fig.add_axes(ax)533 p1.show()534 surfacefig = 'solution_surface%s' % leg_label[j]535 p1.savefig(surfacefig)536 p1.close()537 538 #### finished generating quantities for all swwfiles #####539 540 # x profile for given time541 if surface is True:542 figure(11)543 plot(eastings[tindex,:,j], stages[tindex,:,j])544 xlabel('x')545 ylabel('stage')546 profilefig = 'solution_xprofile'547 savefig('profilefig')548 549 elev_output = []550 if generate_fig is True:551 depth_axis = axis([starttime/scale, time_max/scale, -0.1,552 max(max_depths)*1.1])553 stage_axis = axis([starttime/scale, time_max/scale,554 min(min_stages), max(max_stages)*1.1])555 vel_axis = axis([starttime/scale, time_max/scale,556 min(min_speeds), max(max_speeds)*1.1])557 mom_axis = axis([starttime/scale, time_max/scale,558 min(min_momentums), max(max_momentums)*1.1])559 xmom_axis = axis([starttime/scale, time_max/scale,560 min(min_xmomentums), max(max_xmomentums)*1.1])561 ymom_axis = axis([starttime/scale, time_max/scale,562 min(min_ymomentums), max(max_ymomentums)*1.1])563 cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']564 nn = len(plot_quantity)565 no_cols = 2566 567 if len(label_id) > 1: graphname_report = []568 pp = 1569 div = 11.570 cc = 0571 for k in gauge_index:572 g = gauges[k]573 count1 = 0574 if report == True and len(label_id) > 1:575 s = '\\begin{figure}[ht] \n' \576 '\\centering \n' \577 '\\begin{tabular}{cc} \n'578 fid.write(s)579 if len(label_id) > 1: graphname_report = []580 581 #### generate figures for each gauge ####582 for j, f in enumerate(f_list):583 ion()584 hold(True)585 count = 0586 where1 = 0587 where2 = 0588 word_quantity = ''589 if report == True and len(label_id) == 1:590 s = '\\begin{figure}[hbt] \n' \591 '\\centering \n' \592 '\\begin{tabular}{cc} \n'593 fid.write(s)594 595 for which_quantity in plot_quantity:596 count += 1597 where1 += 1598 figure(count, frameon = False)599 if which_quantity == 'depth':600 plot(model_time[0:n[j]-1,k,j],601 depths[0:n[j]-1,k,j], '-', c = cstr[j])602 units = 'm'603 axis(depth_axis)604 if which_quantity == 'stage':605 if elevations[0,k,j] <= 0:606 plot(model_time[0:n[j]-1,k,j],607 stages[0:n[j]-1,k,j], '-', c = cstr[j])608 axis(stage_axis)609 else:610 plot(model_time[0:n[j]-1,k,j],611 depths[0:n[j]-1,k,j], '-', c = cstr[j])612 #axis(depth_axis)613 units = 'm'614 if which_quantity == 'momentum':615 plot(model_time[0:n[j]-1,k,j],616 momenta[0:n[j]-1,k,j], '-', c = cstr[j])617 axis(mom_axis)618 units = 'm^2 / sec'619 if which_quantity == 'xmomentum':620 plot(model_time[0:n[j]-1,k,j],621 xmom[0:n[j]-1,k,j], '-', c = cstr[j])622 axis(xmom_axis)623 units = 'm^2 / sec'624 if which_quantity == 'ymomentum':625 plot(model_time[0:n[j]-1,k,j],626 ymom[0:n[j]-1,k,j], '-', c = cstr[j])627 axis(ymom_axis)628 units = 'm^2 / sec'629 if which_quantity == 'speed':630 plot(model_time[0:n[j]-1,k,j],631 speed[0:n[j]-1,k,j], '-', c = cstr[j])632 axis(vel_axis)633 units = 'm / sec'634 if which_quantity == 'bearing':635 plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',636 model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.',637 model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')638 units = 'degrees from North'639 #ax = axis([time_min, time_max, 0.0, 360.0])640 legend(('Bearing','West','East'))641 642 if time_unit is 'mins': xlabel('time (mins)')643 if time_unit is 'hours': xlabel('time (hours)')644 #if which_quantity == 'stage' \645 # and elevations[0:n[j]-1,k,j] > 0:646 # ylabel('%s (%s)' %('depth', units))647 #else:648 # ylabel('%s (%s)' %(which_quantity, units))649 #ylabel('%s (%s)' %('wave height', units))650 ylabel('%s (%s)' %(which_quantity, units))651 if len(label_id) > 1: legend((leg_label),loc='upper right')652 653 #gaugeloc1 = gaugeloc.replace(' ','')654 #gaugeloc2 = gaugeloc1.replace('_','')655 gaugeloc2 = str(locations[k]).replace(' ','')656 graphname = '%sgauge%s_%s' %(file_loc[j],657 gaugeloc2,658 which_quantity)659 660 if report == True and len(label_id) > 1:661 figdir = getcwd()+sep+'report_figures'+sep662 if access(figdir,F_OK) == 0 :663 mkdir (figdir)664 latex_file_loc = figdir.replace(sep,altsep)665 # storing files in production directory666 graphname_latex = '%sgauge%s%s' \667 % (latex_file_loc, gaugeloc2,668 which_quantity)669 # giving location in latex output file670 graphname_report_input = '%sgauge%s%s' % \671 ('..' + altsep +672 'report_figures' + altsep,673 gaugeloc2, which_quantity)674 graphname_report.append(graphname_report_input)675 676 # save figures in production directory for report677 savefig(graphname_latex)678 679 if report == True:680 figdir = getcwd() + sep + 'report_figures' + sep681 if access(figdir,F_OK) == 0:682 mkdir(figdir)683 latex_file_loc = figdir.replace(sep,altsep)684 685 if len(label_id) == 1:686 # storing files in production directory687 graphname_latex = '%sgauge%s%s%s' % \688 (latex_file_loc, gaugeloc2,689 which_quantity, label_id2)690 # giving location in latex output file691 graphname_report = '%sgauge%s%s%s' % \692 ('..' + altsep +693 'report_figures' + altsep,694 gaugeloc2, which_quantity,695 label_id2)696 s = '\includegraphics' \697 '[width=0.49\linewidth, height=50mm]{%s%s}' % \698 (graphname_report, '.png')699 fid.write(s)700 if where1 % 2 == 0:701 s = '\\\\ \n'702 where1 = 0703 else:704 s = '& \n'705 fid.write(s)706 savefig(graphname_latex)707 708 if title_on == True:709 title('%s scenario: %s at %s gauge' % \710 (label_id, which_quantity, gaugeloc2))711 #title('Gauge %s (MOST elevation %.2f, ' \712 # 'ANUGA elevation %.2f)' % \713 # (gaugeloc2, elevations[10,k,0],714 # elevations[10,k,1]))715 716 savefig(graphname) # save figures with sww file717 718 if report == True and len(label_id) == 1:719 for i in range(nn-1):720 if nn > 2:721 if plot_quantity[i] == 'stage' \722 and elevations[0,k,j] > 0:723 word_quantity += 'depth' + ', '724 else:725 word_quantity += plot_quantity[i] + ', '726 else:727 if plot_quantity[i] == 'stage' \728 and elevations[0,k,j] > 0:729 word_quantity += 'depth' + ', '730 else:731 word_quantity += plot_quantity[i]732 733 if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:734 word_quantity += ' and ' + 'depth'735 else:736 word_quantity += ' and ' + plot_quantity[nn-1]737 caption = 'Time series for %s at %s location ' \738 '(elevation %.2fm)' % \739 (word_quantity, locations[k], elev[k])740 if elev[k] == 0.0:741 caption = 'Time series for %s at %s location ' \742 '(elevation %.2fm)' % \743 (word_quantity, locations[k],744 elevations[0,k,j])745 east = gauges[0]746 north = gauges[1]747 elev_output.append([locations[k], east, north,748 elevations[0,k,j]])749 label = '%sgauge%s' % (label_id2, gaugeloc2)750 s = '\end{tabular} \n' \751 '\\caption{%s} \n' \752 '\label{fig:%s} \n' \753 '\end{figure} \n \n' % (caption, label)754 fid.write(s)755 cc += 1756 if cc % 6 == 0: fid.write('\\clearpage \n')757 savefig(graphname_latex)758 759 if report == True and len(label_id) > 1:760 for i in range(nn-1):761 if nn > 2:762 if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:763 word_quantity += 'depth' + ','764 else:765 word_quantity += plot_quantity[i] + ', '766 else:767 if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:768 word_quantity += 'depth'769 else:770 word_quantity += plot_quantity[i]771 where1 = 0772 count1 += 1773 index = j*len(plot_quantity)774 for which_quantity in plot_quantity:775 where1 += 1776 s = '\includegraphics' \777 '[width=0.49\linewidth, height=50mm]{%s%s}' % \778 (graphname_report[index], '.png')779 index += 1780 fid.write(s)781 if where1 % 2 == 0:782 s = '\\\\ \n'783 where1 = 0784 else:785 s = '& \n'786 fid.write(s)787 word_quantity += ' and ' + plot_quantity[nn-1]788 label = 'gauge%s' %(gaugeloc2)789 caption = 'Time series for %s at %s location ' \790 '(elevation %.2fm)' % \791 (word_quantity, locations[k], elev[k])792 if elev[k] == 0.0:793 caption = 'Time series for %s at %s location ' \794 '(elevation %.2fm)' % \795 (word_quantity, locations[k],796 elevations[0,k,j])797 thisgauge = gauges[k]798 east = thisgauge[0]799 north = thisgauge[1]800 elev_output.append([locations[k], east, north,801 elevations[0,k,j]])802 803 s = '\end{tabular} \n' \804 '\\caption{%s} \n' \805 '\label{fig:%s} \n' \806 '\end{figure} \n \n' % (caption, label)807 fid.write(s)808 if float((k+1)/div - pp) == 0.:809 fid.write('\\clearpage \n')810 pp += 1811 #### finished generating figures ###812 813 close('all')814 815 return texfile2, elev_output816 319 817 320
Note: See TracChangeset
for help on using the changeset viewer.
