Changeset 4876


Ignore:
Timestamp:
Dec 6, 2007, 5:37:13 PM (17 years ago)
Author:
nick
Message:

added csv2timeseries_graphs. This creates graphs from csv files. Can be used in combination with gauges_sww2csv.

Location:
anuga_core/source/anuga/abstract_2d_finite_volumes
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

    r4647 r4876  
    1515
    1616from sys import platform
     17
     18from anuga.pmesh.mesh import Mesh
     19from anuga.shallow_water import Domain, Transmissive_boundary
     20from anuga.shallow_water.data_manager import get_dataobject
     21from csv import reader,writer
     22import time
     23import string
    1724
    1825def test_function(x, y):
     
    13921399       
    13931400
    1394    
     1401    def test_gauges_sww2csv(self):
     1402
     1403        def elevation_function(x, y):
     1404            return -x
     1405       
     1406        """Most of this test was copied from test_interpolate test_interpole_sww2csv
     1407       
     1408        This is testing the gauge_sww2csv function, by creating a sww file and
     1409        then exporting the gauges and checking the results.
     1410        """
     1411       
     1412        # create mesh
     1413        mesh_file = tempfile.mktemp(".tsh")   
     1414        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
     1415        m = Mesh()
     1416        m.add_vertices(points)
     1417        m.auto_segment()
     1418        m.generate_mesh(verbose=False)
     1419        m.export_mesh_file(mesh_file)
     1420       
     1421        #Create shallow water domain
     1422        domain = Domain(mesh_file)
     1423        os.remove(mesh_file)
     1424       
     1425        domain.default_order=2
     1426        domain.beta_h = 0
     1427
     1428        #Set some field values
     1429        domain.set_quantity('elevation', elevation_function)
     1430        domain.set_quantity('friction', 0.03)
     1431        domain.set_quantity('xmomentum', 3.0)
     1432        domain.set_quantity('ymomentum', 4.0)
     1433
     1434        ######################
     1435        # Boundary conditions
     1436        B = Transmissive_boundary(domain)
     1437        domain.set_boundary( {'exterior': B})
     1438
     1439        # This call mangles the stage values.
     1440        domain.distribute_to_vertices_and_edges()
     1441        domain.set_quantity('stage', 1.0)
     1442
     1443
     1444        domain.set_name('datatest' + str(time.time()))
     1445        domain.format = 'sww'
     1446        domain.smooth = True
     1447        domain.reduction = mean
     1448
     1449        sww = get_dataobject(domain)
     1450        sww.store_connectivity()
     1451        sww.store_timestep(['stage', 'xmomentum', 'ymomentum','elevation'])
     1452        domain.set_quantity('stage', 10.0) # This is automatically limited
     1453        # so it will not be less than the elevation
     1454        domain.time = 2.
     1455        sww.store_timestep(['stage','elevation', 'xmomentum', 'ymomentum'])
     1456
     1457        # test the function
     1458        points = [[5.0,1.],[0.5,2.]]
     1459
     1460        points_file = tempfile.mktemp(".csv")
     1461#        points_file = 'test_point.csv'
     1462        file_id = open(points_file,"w")
     1463        file_id.write("name, easting, northing, elevation \n\
     1464point1, 5.0, 1.0, 3.0\n\
     1465point2, 0.5, 2.0, 9.0\n")
     1466        file_id.close()
     1467
     1468       
     1469        gauges_sww2csv(sww.filename,
     1470                            points_file,
     1471                            verbose=False,
     1472                            use_cache=False)
     1473
     1474        point1_answers_array = [[0.0,1.0,-5.0,3.0,4.0], [2.0,10.0,-5.0,3.0,4.0]]
     1475        point1_filename = 'gauge_point1.csv'
     1476        point1_handle = file(point1_filename)
     1477        point1_reader = reader(point1_handle)
     1478        point1_reader.next()
     1479
     1480        line=[]
     1481        for i,row in enumerate(point1_reader):
     1482#            print 'i',i,'row',row
     1483            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4])])
     1484            #print 'assert line',line[i],'point1',point1_answers_array[i]
     1485            assert allclose(line[i], point1_answers_array[i])
     1486
     1487        point2_answers_array = [[0.0,1.0,-0.5,3.0,4.0], [2.0,10.0,-0.5,3.0,4.0]]
     1488        point2_filename = 'gauge_point2.csv'
     1489        point2_handle = file(point2_filename)
     1490        point2_reader = reader(point2_handle)
     1491        point2_reader.next()
     1492                       
     1493        line=[]
     1494        for i,row in enumerate(point2_reader):
     1495#            print 'i',i,'row',row
     1496            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4])])
     1497            #print 'assert line',line[i],'point1',point1_answers_array[i]
     1498            assert allclose(line[i], point2_answers_array[i])
     1499                         
     1500        # clean up
     1501        point1_handle.close()
     1502        point2_handle.close()
     1503        #print "sww.filename",sww.filename
     1504        os.remove(sww.filename)
     1505        os.remove(points_file)
     1506        os.remove(point1_filename)
     1507        os.remove(point2_filename)
     1508
     1509
     1510
     1511    def test_gauges_sww2csv1(self):
     1512        from anuga.pmesh.mesh import Mesh
     1513        from anuga.shallow_water import Domain, Transmissive_boundary
     1514        from anuga.shallow_water.data_manager import get_dataobject
     1515        from csv import reader,writer
     1516        import time
     1517        import string
     1518
     1519        def elevation_function(x, y):
     1520            return -x
     1521       
     1522        """Most of this test was copied from test_interpolate test_interpole_sww2csv
     1523       
     1524        This is testing the gauge_sww2csv function, by creating a sww file and
     1525        then exporting the gauges and checking the results.
     1526       
     1527        This tests the ablity not to have elevation in the points file and
     1528        not store xmomentum and ymomentum
     1529        """
     1530       
     1531        # create mesh
     1532        mesh_file = tempfile.mktemp(".tsh")   
     1533        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
     1534        m = Mesh()
     1535        m.add_vertices(points)
     1536        m.auto_segment()
     1537        m.generate_mesh(verbose=False)
     1538        m.export_mesh_file(mesh_file)
     1539       
     1540        #Create shallow water domain
     1541        domain = Domain(mesh_file)
     1542        os.remove(mesh_file)
     1543       
     1544        domain.default_order=2
     1545        domain.beta_h = 0
     1546
     1547        #Set some field values
     1548        domain.set_quantity('elevation', elevation_function)
     1549        domain.set_quantity('friction', 0.03)
     1550        domain.set_quantity('xmomentum', 3.0)
     1551        domain.set_quantity('ymomentum', 4.0)
     1552
     1553        ######################
     1554        # Boundary conditions
     1555        B = Transmissive_boundary(domain)
     1556        domain.set_boundary( {'exterior': B})
     1557
     1558        # This call mangles the stage values.
     1559        domain.distribute_to_vertices_and_edges()
     1560        domain.set_quantity('stage', 1.0)
     1561
     1562
     1563        domain.set_name('datatest' + str(time.time()))
     1564        domain.format = 'sww'
     1565        domain.smooth = True
     1566        domain.reduction = mean
     1567
     1568        sww = get_dataobject(domain)
     1569        sww.store_connectivity()
     1570        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     1571        domain.set_quantity('stage', 10.0) # This is automatically limited
     1572        # so it will not be less than the elevation
     1573        domain.time = 2.
     1574        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     1575
     1576        # test the function
     1577        points = [[5.0,1.],[0.5,2.]]
     1578
     1579        points_file = tempfile.mktemp(".csv")
     1580#        points_file = 'test_point.csv'
     1581        file_id = open(points_file,"w")
     1582        file_id.write("name, easting, northing \n\
     1583point1, 5.0, 1.0\n\
     1584point2, 0.5, 2.0\n")
     1585        file_id.close()
     1586
     1587        gauges_sww2csv(sww.filename,
     1588                            points_file,
     1589                            quantities=['Stage', 'elevation'],
     1590                            use_cache=False,
     1591                            verbose=False)
     1592
     1593        point1_answers_array = [[0.0,1.0,-5.0], [2.0,10.0,-5.0]]
     1594        point1_filename = 'gauge_point1.csv'
     1595        point1_handle = file(point1_filename)
     1596        point1_reader = reader(point1_handle)
     1597        point1_reader.next()
     1598
     1599        line=[]
     1600        for i,row in enumerate(point1_reader):
     1601#            print 'i',i,'row',row
     1602            line.append([float(row[0]),float(row[1]),float(row[2])])
     1603            #print 'line',line[i],'point1',point1_answers_array[i]
     1604            assert allclose(line[i], point1_answers_array[i])
     1605
     1606        point2_answers_array = [[0.0,1.0,-0.5], [2.0,10.0,-0.5]]
     1607        point2_filename = 'gauge_point2.csv'
     1608        point2_handle = file(point2_filename)
     1609        point2_reader = reader(point2_handle)
     1610        point2_reader.next()
     1611                       
     1612        line=[]
     1613        for i,row in enumerate(point2_reader):
     1614#            print 'i',i,'row',row
     1615            line.append([float(row[0]),float(row[1]),float(row[2])])
     1616#            print 'line',line[i],'point1',point1_answers_array[i]
     1617            assert allclose(line[i], point2_answers_array[i])
     1618                         
     1619        # clean up
     1620        point1_handle.close()
     1621        point2_handle.close()
     1622        #print "sww.filename",sww.filename
     1623        os.remove(sww.filename)
     1624        os.remove(points_file)
     1625        os.remove(point1_filename)
     1626        os.remove(point2_filename)
    13951627       
    13961628#-------------------------------------------------------------
    13971629if __name__ == "__main__":
    13981630    suite = unittest.makeSuite(Test_Util,'test')
    1399 #    suite = unittest.makeSuite(Test_Util,'test_get_min_max_values')
    1400     runner = unittest.TextTestRunner(verbosity=2)
     1631#    suite = unittest.makeSuite(Test_Util,'test_gauges_sww')
     1632#    runner = unittest.TextTestRunner(verbosity=2)
     1633    runner = unittest.TextTestRunner(verbosity=1)
    14011634    runner.run(suite)
    14021635
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r4820 r4876  
    99import os
    1010
    11 from os import remove, mkdir, access, F_OK, W_OK, sep,mkdir
     11from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep,mkdir
    1212from os.path import exists, basename, split,join
    1313from warnings import warn
     
    14941494    return xc
    14951495
    1496 def make_plots_from_csv_file(directories_dic={},
     1496def csv2timeseries_graphs(directories_dic={},
    14971497                            output_dir='',
    14981498                            base_name=None,
    14991499                            plot_numbers='0',
    1500                             quantities=['Stage'],
     1500                            quantities=['stage'],
    15011501                            extra_plot_name='',
    15021502                            assess_all_csv_files=True,                           
     
    15041504                            verbose=False):
    15051505                               
    1506     """Read in csv files that have the right header information and
     1506    """
     1507    WARNING!! NO UNIT TESTS... could make some as to test filename..but have
     1508    spend ages on this already...
     1509   
     1510    Read in csv files that have the right header information and
    15071511    plot time series such as Stage, Speed, etc. Will also plot several
    15081512    time series on one plot. Filenames must follow this convention,
     
    15651569    #be very good and very flexable.... but this works for now!
    15661570
    1567     #FIXME plot all available data from the csv file, currently will
    1568     #only plot Stage,Speed and Momentum and only if the header is
    1569     #named like wise
    1570 
    1571 #    quantities_dic={'time':'time (hour)'}
    1572 #    if 'time' in quantities:
    1573 #        quantities_dic['time']='time (hour)'
    1574 #    if 'Time' in quantities:
    1575 #        quantities_dic['Time']='time (hour)'
    1576 #    if 'stage' in quantities:
    1577 #        quantities_dic['stage']='wave height (m)'
    1578 #    if 'Stage' in quantities:
    1579 #        quantities_dic['Stage']='wave height (m)'
    1580 #    if 'speed' in quantities:
    1581 #        quantities_dic['speed']='speed (m/s)'
    1582 #    if 'Speed' in quantities:
    1583 #        quantities_dic['Speed']='speed (m/s)'
    1584 #    if 'stage' or 'Stage' in quantities:
    1585 #        quantities_dic['stage']='speed (m/s)'
    1586 #    if 'stage' or 'Stage' in quantities:
    1587 #        quantities_dic['stage']='speed (m/s)'
     1571    #FIXME MAKE NICER at the end i have started a new version...
     1572    # it almost completely works...
    15881573
    15891574    seconds_in_hour = 3600
     
    15921577    speed_label = 'speed (m/s)'
    15931578    momentum_label = 'momentum (m^2/sec)'
     1579    xmomentum_label = 'momentum (m^2/sec)'
     1580    ymomentum_label = 'momentum (m^2/sec)'
     1581    depth_label = 'water depth (m)'
     1582    bearing_label = 'degrees (o)'
    15941583   
    15951584    if extra_plot_name != '':
     
    16061595
    16071596    #use all the files to get the values for the plot axis
    1608     max_st=max_sp=max_mom=min_st=min_sp=min_mom=max_t=min_t=0.
     1597    max_xmom=min_xmom=max_ymom=min_ymom=max_st=max_sp=max_mom=min_st=min_sp=min_mom=\
     1598    max_depth=min_depth=max_bearing=min_bearing=max_t=min_t=0.
    16091599    max_start_time= 0.
    16101600    min_start_time = 100000
     
    16471637            directory_add_tide = directories_dic[directory][2]
    16481638
    1649             time = [float(x) for x in attribute_dic["Time"]]
     1639            time = [float(x) for x in attribute_dic["time"]]
    16501640            min_t, max_t = get_min_max_values(time,min_t,max_t)
    16511641           
    1652             stage = [float(x) for x in attribute_dic["Stage"]]
     1642            stage = [float(x) for x in attribute_dic["stage"]]
    16531643            stage =array(stage)+directory_add_tide
    16541644            min_st, max_st = get_min_max_values(stage,min_st,max_st)
     1645       
     1646            if "speed" in quantities:
     1647                speed = [float(x) for x in attribute_dic["speed"]]
     1648                min_sp, max_sp = get_min_max_values(speed,min_sp,max_sp)
    16551649           
    1656             speed = [float(x) for x in attribute_dic["Speed"]]
    1657             min_sp, max_sp = get_min_max_values(speed,min_sp,max_sp)
    1658            
    1659             momentum = [float(x) for x in attribute_dic["Momentum"]]
    1660             min_mom, max_mom = get_min_max_values(momentum,
     1650            if "momentum" in quantities:
     1651                momentum = [float(x) for x in attribute_dic["momentum"]]
     1652                min_mom, max_mom = get_min_max_values(momentum,
    16611653                                                  min_mom,
    16621654                                                  max_mom)
     1655            if "xmomentum" in quantities:
     1656                xmomentum = [float(x) for x in attribute_dic["xmomentum"]]
     1657                min_xmom, max_xmom = get_min_max_values(xmomentum,
     1658                                                  min_xmom,
     1659                                                  max_xmom)
     1660            if "ymomentum" in quantities:
     1661                ymomentum = [float(x) for x in attribute_dic["ymomentum"]]
     1662                min_ymom, max_ymom = get_min_max_values(ymomentum,
     1663                                                  min_ymom,
     1664                                                  max_ymom)
     1665            if "depth" in quantities:
     1666                depth = [float(x) for x in attribute_dic["depth"]]
     1667                min_depth, max_depth = get_min_max_values(depth,
     1668                                                  min_depth,
     1669                                                  max_depth)
     1670            if "bearing" in quantities:
     1671                bearing = [float(x) for x in attribute_dic["bearing"]]
     1672                min_bearing, max_bearing = get_min_max_values(bearing,
     1673                                                  min_bearing,
     1674                                                  max_bearing)
    16631675                                                                     
    16641676#            print 'min_sp, max_sp',min_sp, max_sp,
     
    16761688                 (max_t+max_start_time)/seconds_in_hour,
    16771689                 min_sp, max_sp)
     1690    xmomentum_axis = (min_start_time/seconds_in_hour,
     1691                    (max_t+max_start_time)/seconds_in_hour,
     1692                     min_xmom, max_xmom)
     1693    ymomentum_axis = (min_start_time/seconds_in_hour,
     1694                    (max_t+max_start_time)/seconds_in_hour,
     1695                     min_ymom, max_ymom)
    16781696    momentum_axis = (min_start_time/seconds_in_hour,
    16791697                    (max_t+max_start_time)/seconds_in_hour,
    16801698                     min_mom, max_mom)
    1681 #    ion()
    1682 #    pylab.hold()
    1683    
    1684    
     1699    depth_axis = (min_start_time/seconds_in_hour,
     1700                    (max_t+max_start_time)/seconds_in_hour,
     1701                     min_depth, max_depth)
     1702    bearing_axis = (min_start_time/seconds_in_hour,
     1703                    (max_t+max_start_time)/seconds_in_hour,
     1704                     min_bearing, max_bearing)
     1705 
    16851706    cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k']
    16861707   
     
    17041725            attribute_dic, title_index_dic = csv2dict(file+'.csv')
    17051726            #get data from dict in to list
    1706             t = [float(x) for x in attribute_dic["Time"]]
     1727            t = [float(x) for x in attribute_dic["time"]]
    17071728
    17081729            #do maths to list by changing to array
    17091730            t=(array(t)+directory_start_time)/seconds_in_hour
    1710             stage = [float(x) for x in attribute_dic["Stage"]]
    1711             speed = [float(x) for x in attribute_dic["Speed"]]
    1712             momentum = [float(x) for x in attribute_dic["Momentum"]]
     1731            stage = [float(x) for x in attribute_dic["stage"]]
     1732            if "speed" in quantities:
     1733                speed = [float(x) for x in attribute_dic["speed"]]
     1734            if "xmomentum" in quantities:
     1735                xmomentum = [float(x) for x in attribute_dic["xmomentum"]]
     1736            if "ymomentum" in quantities:
     1737                ymomentum = [float(x) for x in attribute_dic["ymomentum"]]
     1738            if "momentum" in quantities:
     1739                momentum = [float(x) for x in attribute_dic["momentum"]]
     1740            if "depth" in quantities:
     1741                depth = [float(x) for x in attribute_dic["depth"]]
     1742            if "bearing" in quantities:
     1743                bearing = [float(x) for x in attribute_dic["bearing"]]
    17131744                       
    17141745            #Can add tide so plots are all on the same tide height
     
    17181749            max_ele=-100000
    17191750            min_ele=100000
    1720             elevation = [float(x) for x in attribute_dic["Elevation"]]
     1751            elevation = [float(x) for x in attribute_dic["elevation"]]
    17211752           
    17221753            min_ele, max_ele = get_min_max_values(elevation,
     
    17431774                    #print k,l, legend_list_dic[j]
    17441775         
    1745             if "Stage" in quantities:
     1776            if "stage" in quantities:
    17461777                pylab.figure(100+j)
    17471778                pylab.plot(t, stage, c = cstr[i], linewidth=1)
     
    17501781                pylab.axis(stage_axis)
    17511782                pylab.legend(legend_list,loc='upper right')
    1752                 figname = '%sstage_%s%s.png' %(output_dir+sep,
     1783                if output_dir =='':
     1784                    figname = '%sstage_%s%s.png' %(directory+sep,
    17531785                                               base_name+number,
    17541786                                               extra_plot_name)
     1787                else:   
     1788                    figname = '%sstage_%s%s.png' %(output_dir+sep,
     1789                                               base_name+number,
     1790                                               extra_plot_name)
     1791                if verbose: print 'saving figure here %s' %figname
    17551792                pylab.savefig(figname)
    1756             if "Speed" in quantities:
     1793
     1794            if "speed" in quantities:
    17571795                pylab.figure(200+j)
    17581796                pylab.plot(t, speed, c = cstr[i], linewidth=1)
     
    17611799                pylab.axis(speed_axis)
    17621800                pylab.legend(legend_list,loc='upper right')
    1763                 figname = '%sspeed_%s%s.png' %(output_dir+sep,
     1801                if output_dir =='':
     1802                    figname = '%sspeed_%s%s.png' %(directory+sep,
    17641803                                               base_name+number,
    17651804                                               extra_plot_name)
     1805                else:
     1806                    figname = '%sspeed_%s%s.png' %(output_dir+sep,
     1807                                               base_name+number,
     1808                                               extra_plot_name)
     1809                if verbose: print 'saving figure here %s' %figname
    17661810                pylab.savefig(figname)
    1767             if "Momentum" in quantities:
     1811            if "xmomentum" in quantities:
    17681812                pylab.figure(300+j)
     1813                pylab.plot(t, xmomentum, c = cstr[i], linewidth=1)
     1814                pylab.xlabel(time_label)
     1815                pylab.ylabel(xmomentum_label)
     1816                pylab.axis(xmomentum_axis)
     1817                pylab.legend(legend_list,loc='upper right')
     1818                if output_dir =='':
     1819                    figname = '%sxmomentum_%s%s.png' %(directory+sep,
     1820                                                  base_name+number,
     1821                                                  extra_plot_name)
     1822                else:
     1823                    figname = '%sxmomentum_%s%s.png' %(output_dir+sep,
     1824                                                  base_name+number,
     1825                                                  extra_plot_name)
     1826                if verbose: print 'saving figure here %s' %figname
     1827                pylab.savefig(figname)
     1828               
     1829            if "ymomentum" in quantities:
     1830                pylab.figure(400+j)
     1831                pylab.plot(t, ymomentum, c = cstr[i], linewidth=1)
     1832                pylab.xlabel(time_label)
     1833                pylab.ylabel(ymomentum_label)
     1834                pylab.axis(ymomentum_axis)
     1835                pylab.legend(legend_list,loc='upper right')
     1836                if output_dir =='':
     1837                    figname = '%symomentum_%s%s.png' %(directory+sep,
     1838                                                  base_name+number,
     1839                                                  extra_plot_name)
     1840                else:
     1841                    figname = '%symomentum_%s%s.png' %(output_dir+sep,
     1842                                                  base_name+number,
     1843                                                  extra_plot_name)
     1844                if verbose: print 'saving figure here %s' %figname
     1845                pylab.savefig(figname)
     1846
     1847            if "momentum" in quantities:
     1848                pylab.figure(500+j)
    17691849                pylab.plot(t, momentum, c = cstr[i], linewidth=1)
    17701850                pylab.xlabel(time_label)
     
    17721852                pylab.axis(momentum_axis)
    17731853                pylab.legend(legend_list,loc='upper right')
    1774                 figname = '%smomentum_%s%s.png' %(output_dir+sep,
     1854                if output_dir =='':
     1855                    figname = '%smomentum_%s%s.png' %(directory+sep,
    17751856                                                  base_name+number,
    17761857                                                  extra_plot_name)
     1858                else:
     1859                    figname = '%smomentum_%s%s.png' %(output_dir+sep,
     1860                                                  base_name+number,
     1861                                                  extra_plot_name)
     1862                if verbose: print 'saving figure here %s' %figname
    17771863                pylab.savefig(figname)
     1864
     1865            if "bearing" in quantities:
     1866                pylab.figure(600+j)
     1867                pylab.plot(t, bearing, c = cstr[i], linewidth=1)
     1868                pylab.xlabel(time_label)
     1869                pylab.ylabel(bearing_label)
     1870                pylab.axis(bearing_axis)
     1871                pylab.legend(legend_list,loc='upper right')
     1872                if output_dir =='':
     1873                    figname = '%sbearing_%s%s.png' %(output_dir+sep,
     1874                                                  base_name+number,
     1875                                                  extra_plot_name)
     1876                else:
     1877                    figname = '%sbearing_%s%s.png' %(output_dir+sep,
     1878                                                  base_name+number,
     1879                                                  extra_plot_name)
     1880                if verbose: print 'saving figure here %s' %figname
     1881                pylab.savefig(figname)
     1882
     1883            if "depth" in quantities:
     1884                pylab.figure(700+j)
     1885                pylab.plot(t, depth, c = cstr[i], linewidth=1)
     1886                pylab.xlabel(time_label)
     1887                pylab.ylabel(depth_label)
     1888                pylab.axis(depth_axis)
     1889                pylab.legend(legend_list,loc='upper right')
     1890                if output_dir == '':
     1891                    figname = '%sdepth_%s%s.png' %(directory,
     1892                                                  base_name+number,
     1893                                                  extra_plot_name)
     1894                else:
     1895                    figname = '%sdepth_%s%s.png' %(output_dir+sep,
     1896                                                  base_name+number,
     1897                                                  extra_plot_name)
     1898                if verbose: print 'saving figure here %s' %figname
     1899                pylab.savefig(figname)
     1900               
    17781901    if verbose: print 'Closing all plots'
    17791902    pylab.close('all')
    17801903    del pylab
    17811904    if verbose: print 'Finished closing plots'
     1905   
     1906def XXX_csv2timeseries_graphs(directories_dic={},
     1907                            output_dir='',
     1908                            base_name='',
     1909                            plot_numbers='',
     1910                            quantities=['stage'],
     1911                            extra_plot_name='',
     1912                            assess_all_csv_files=False,                           
     1913                            create_latex=False,
     1914                            verbose=False):
     1915                               
     1916    """     
     1917    NOTE! DOES not work 100% yet. eg don't use assess_all_csv_files=True
     1918    everything else is ok?                 
     1919    Use sww2timeseries....           
     1920   
     1921    Read in csv files that have the right header information and
     1922    plot time series such as Stage, Speed, etc. Will also plot several
     1923    time series on one plot. Filenames must follow this convention,
     1924    <base_name><plot_number>.csv eg gauge_timeseries3.csv
     1925
     1926    Each file represents a location and within each file there are
     1927    time, quantity columns.
     1928   
     1929    For example:   
     1930    if "directories_dic" defines 4 directories and in each directories
     1931    there is a csv files corresponding to the right "plot_numbers",
     1932    this will create a plot with 4 lines one for each directory AND
     1933    one plot for each "quantities".  ??? FIXME: unclear.
     1934   
     1935    Inputs:
     1936        directories_dic: dictionary of directory with values (plot
     1937                         legend name for directory), (start time of
     1938                         the time series) and the (value to add to
     1939                         stage if needed). For example
     1940                        {dir1:['Anuga_ons',5000, 0],
     1941                         dir2:['b_emoth',5000,1.5],
     1942                         dir3:['b_ons',5000,1.5]}
     1943                         
     1944        output_dir: directory for the plot outputs
     1945       
     1946        base_name: common name to all the csv files to be read, Note is
     1947                   ignored if assess_all_csv_files=True
     1948       
     1949        plot_numbers: a String list of numbers to plot. For example
     1950                      [0-4,10,15-17] will read and attempt to plot
     1951                       the follow 0,1,2,3,4,10,15,16,17
     1952                       NOTE: if no plot numbers this will create only
     1953                       one plot per quantity
     1954        quantities: Currently limited to "stage", "speed", and
     1955                     "Momentum", should be changed to incorporate
     1956                    any quantity read from CSV file....
     1957                   
     1958        extra_plot_name: A string that is appended to the end of the
     1959                         output filename.
     1960                   
     1961        assess_all_csv_files: if true it will read ALL csv file with
     1962                             "base_name", regardless of 'plot_numbers'
     1963                              and determine a uniform set of axes for
     1964                              Stage, Speed and Momentum. IF FALSE it
     1965                              will only read the csv file within the
     1966                             'plot_numbers'
     1967                             
     1968        create_latex: NOT IMPLEMENTED YET!! sorry Jane....
     1969       
     1970    OUTPUTS: None, it saves the plots to
     1971              <output_dir><base_name><plot_number><extra_plot_name>.png
     1972    """
     1973
     1974    import pylab# import plot, xlabel, ylabel, savefig, \
     1975                #ion, hold, axis, close, figure, legend
     1976    from os import sep
     1977    from anuga.shallow_water.data_manager import \
     1978                               get_all_files_with_extension, csv2dict
     1979    #find all the files that meet the specs
     1980
     1981    seconds_in_hour = 3600
     1982   
     1983    quantities_label={}
     1984    quantities_label['time'] = 'time (hour)'
     1985    quantities_label['stage'] = 'wave height (m)'
     1986    quantities_label['speed'] = 'speed (m/s)'
     1987    quantities_label['momentum'] = 'momentum (m^2/sec)'
     1988    quantities_label['depth'] = 'momentum (m^2/sec)'
     1989    quantities_label['xmomentum'] = 'momentum (m^2/sec)'
     1990    quantities_label['ymomentum'] = 'momentum (m^2/sec)'
     1991    quantities_label['bearing'] = 'degrees'
     1992   
     1993    if extra_plot_name != '':
     1994        extra_plot_name='_'+extra_plot_name
     1995
     1996    new_plot_numbers=[]
     1997    #change plot_numbers to list, eg ['0-4','10']
     1998    #to ['0','1','2','3','4','10']
     1999    for i, num_string in enumerate(plot_numbers):
     2000        if '-' in num_string:
     2001            start = int(num_string[:num_string.rfind('-')])
     2002            end = int(num_string[num_string.rfind('-')+1:])+1
     2003            for x in range(start, end):
     2004                new_plot_numbers.append(str(x))
     2005        else:
     2006            new_plot_numbers.append(num_string)
     2007
     2008    #finds all the files that fit the specs and return a list of them
     2009    #so to help find a uniform max and min for the plots...
     2010    list_filenames=[]
     2011    filenames=[]
     2012    all_csv_filenames=[]
     2013    if verbose: print 'Determining files to access for axes ranges \n'
     2014#    print directories_dic.keys, base_name
     2015 
     2016    for i,directory in enumerate(directories_dic.keys()):
     2017        all_csv_filenames.append(get_all_files_with_extension(directory,
     2018                                  base_name,'.csv'))
     2019        if assess_all_csv_files==True:
     2020            list_filenames.append(get_all_files_with_extension(directory,
     2021                                  base_name,'.csv'))
     2022        else:
     2023            if plot_numbers == '':
     2024                list_filenames.append(base_name)
     2025            else:
     2026                for number in new_plot_numbers:
     2027                    filenames.append(base_name+number)
     2028                list_filenames.append(filenames)
     2029    print "list_filenames", list_filenames
     2030
     2031    #use all the files to get the values for the plot axis
     2032    #max_st=max_sp=max_mom=min_st=min_sp=min_mom=max_t=min_t=0.
     2033    max_t=-100
     2034    min_t=100000
     2035    max_start_time= -1000.
     2036    min_start_time = 100000
     2037   
     2038    if verbose: print 'Determining uniform axes \n'
     2039    #this entire loop is to determine the min and max range for the
     2040    #axes of the plots
     2041
     2042    quantities.insert(0,'time') 
     2043   
     2044    for i, directory in enumerate(directories_dic.keys()):
     2045       
     2046        if assess_all_csv_files==False:
     2047            which_csv_to_assess = list_filenames[i]
     2048        else:
     2049            #gets list of filenames for directory "i"
     2050            which_csv_to_assess = all_csv_filenames[i]
     2051
     2052       
     2053        for j, filename in enumerate(which_csv_to_assess):
     2054#        for j, filename in enumerate(list_filenames[i]):
     2055
     2056            dir_filename=join(directory,filename)
     2057            attribute_dic, title_index_dic = csv2dict(dir_filename+
     2058                                                       '.csv')
     2059
     2060            directory_start_time = directories_dic[directory][1]
     2061            directory_add_tide = directories_dic[directory][2]
     2062
     2063            print 'keys',attribute_dic.keys()
     2064            quantity_value={}
     2065            min_quantity_value={}
     2066            max_quantity_value={}
     2067#            print 'Quantities',quantities, attribute_dic
     2068       
     2069            #add time to get values
     2070            for k, quantity in enumerate(quantities):
     2071                quantity_value[quantity] = [float(x) for x in attribute_dic[quantity]]
     2072                min_quantity_value[quantity], max_quantity_value[quantity] = get_min_max_values(quantity_value[quantity],\
     2073                                                                  min_t,max_t)
     2074
     2075                if quantity == 'stage':
     2076                     quantity_value[quantity]=array(quantity_value[quantity])+directory_add_tide
     2077
     2078            #set the time... ???
     2079            if min_start_time > directory_start_time:
     2080                min_start_time = directory_start_time
     2081            if max_start_time < directory_start_time:
     2082                max_start_time = directory_start_time
     2083#            print 'start_time' , max_start_time, min_start_time
     2084   
     2085
     2086    #final step to unifrom axis for the graphs
     2087    quantities_axis={}
     2088    for i, quantity in enumerate(quantities):
     2089        print 'max_t %s max_start_time, seconds_in_hour',max_t, max_start_time, seconds_in_hour
     2090        quantities_axis[quantity] = (min_start_time/seconds_in_hour,
     2091                                        (max_quantity_value['time']+max_start_time)/seconds_in_hour,
     2092                                         min_quantity_value[quantity],
     2093                                         max_quantity_value[quantity])
     2094        print 'axis for quantity %s are %s' %(quantity, quantities_axis[quantity])
     2095
     2096    cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k']
     2097
     2098   
     2099    if verbose: print 'Now start to plot \n'
     2100   
     2101    i_max = len(directories_dic.keys())
     2102    legend_list_dic =[]
     2103    legend_list =[]
     2104    for i, directory in enumerate(directories_dic.keys()):
     2105        if verbose: print'Plotting in %s' %directory, new_plot_numbers
     2106#        for j, number in enumerate(new_plot_numbers):
     2107        for j, filename in enumerate(list_filenames[i]):
     2108           
     2109#            if verbose: print'Starting %s' %base_name+number 
     2110            if verbose: print'Starting %s' %filename 
     2111            directory_name = directories_dic[directory][0]
     2112            directory_start_time = directories_dic[directory][1]
     2113            directory_add_tide = directories_dic[directory][2]
     2114           
     2115            #create an if about the start time and tide hieght
     2116            #if they don't exist
     2117#            file=directory+sep+base_name+number
     2118            #print 'i %s,j %s, number %s, file %s' %(i,j,number,file)
     2119            attribute_dic, title_index_dic = csv2dict(filename+'.csv')
     2120            #get data from dict in to list
     2121            t = [float(x) for x in attribute_dic["time"]]
     2122
     2123            #do maths to list by changing to array
     2124            t=(array(t)+directory_start_time)/seconds_in_hour
     2125
     2126            #finds the maximum elevation, used only as a test
     2127            # and as info in the graphs
     2128            max_ele=-100000
     2129            min_ele=100000
     2130            elevation = [float(x) for x in attribute_dic["elevation"]]
     2131           
     2132            min_ele, max_ele = get_min_max_values(elevation,
     2133                                                  min_ele,
     2134                                                  max_ele)
     2135            if min_ele != max_ele:
     2136                print "Note! Elevation changes in %s" %dir_filename
     2137
     2138            #populates the legend_list_dic with dir_name and the elevation
     2139            if i==0:
     2140                legend_list_dic.append({directory_name:max_ele})
     2141            else:
     2142                print j,max_ele, directory_name, legend_list_dic
     2143                legend_list_dic[j][directory_name]=max_ele
     2144
     2145            # creates a list for the legend after at "legend_dic" has been fully populated
     2146            # only runs on the last iteration for all the gauges(csv) files
     2147            # empties the list before creating it
     2148            if i==i_max-1:
     2149                legend_list=[]
     2150                for k, l in legend_list_dic[j].iteritems():
     2151                    legend_list.append('%s (elevation = %sm)'%(k,l))
     2152                    #print k,l, legend_list_dic[j]
     2153           
     2154            print 'filename',filename, quantities
     2155            #remove time so it is not plotted!
     2156#            quantities.remove('time')
     2157            for k, quantity in enumerate(quantities):
     2158                if quantity != 'time':
     2159                    quantity_value[quantity] = [float(x) for x in attribute_dic[quantity]]
     2160
     2161                    if quantity=='stage':
     2162                        quantity_value[quantity] =array(quantity_value[quantity])+directory_add_tide
     2163       
     2164            #Can add tide so plots are all on the same tide height
     2165#                print'plot stuff',quantities_label['time']
     2166#                print quantities_label[quantity]
     2167#                print quantities_axis[quantity]
     2168                    num=int(k*100+j)
     2169                    pylab.figure(num)
     2170                    pylab.plot(t, quantity_value[quantity], c = cstr[i], linewidth=1)
     2171                    pylab.xlabel(quantities_label['time'])
     2172                    pylab.ylabel(quantities_label[quantity])
     2173                    pylab.axis(quantities_axis[quantity])
     2174                    pylab.legend(legend_list,loc='upper right')
     2175                    if output_dir == '':
     2176                        figname = '%s%s_%s%s.png' %(directory,
     2177                                            quantity,
     2178                                            filename,
     2179                                            extra_plot_name)
     2180                    else:
     2181                        figname = '%s%s_%s%s.png' %(output_dir,
     2182                                            quantity,
     2183                                            filename,
     2184                                            extra_plot_name)
     2185                    if verbose: print 'saving figure here %s' %figname
     2186                    pylab.savefig(figname)
     2187           
     2188    if verbose: print 'Closing all plots'
     2189    pylab.close('all')
     2190    del pylab
     2191    if verbose: print 'Finished closing plots'
     2192
    17822193
    17832194def get_min_max_values(list=None,min1=100,max1=-100):
     
    18482259        file.close()
    18492260
    1850 
     2261def gauges_sww2csv(sww_file,
     2262                   gauge_file,
     2263                   quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
     2264                   verbose=False,
     2265                   use_cache = True):
     2266    """
     2267   
     2268    Inputs:
     2269       
     2270        sww_file: ...
     2271       
     2272        points_file: Assumes that it follows this format
     2273            name, easting, northing, elevation
     2274            point1, 100.3, 50.2, 10.0
     2275            point2, 10.3, 70.3, 78.0
     2276       
     2277    Outputs:
     2278        one file for each gauge/point location in the points file. They
     2279        will be named with this format
     2280            gauge_<name>.csv    eg gauge_point1.csv
    18512281           
    1852 
    1853    
    1854        
     2282        They will all have a header
     2283   
     2284    Interpolate the quantities at a given set of locations, given
     2285    an sww file.
     2286    The results are written to a csv file.
     2287
     2288    In the future let points be a points file.
     2289    And the user choose the quantities.
     2290
     2291    This is currently quite specific.
     2292    If it need to be more general, chagne things.
     2293
     2294    This is really returning speed, not velocity.
     2295    """
     2296    from csv import reader,writer
     2297    from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN
     2298    from Numeric import array, resize,shape,Float,zeros,take,argsort,argmin
     2299    import string
     2300    from anuga.shallow_water.data_manager import get_all_swwfiles
     2301#    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
     2302    #print "points",points
     2303
     2304    assert type(gauge_file) == type(''),\
     2305           'Gauge filename must be a string'
     2306   
     2307    try:
     2308    #    fid = open(gauge_file)
     2309        point_reader = reader(file(gauge_file))
     2310
     2311    except Exception, e:
     2312        msg = 'File "%s" could not be opened: Error="%s"'\
     2313                  %(gauge_file, e)
     2314        raise msg
     2315    if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename
     2316   
     2317   
     2318    point_reader = reader(file(gauge_file))
     2319    points = []
     2320    point_name = []
     2321   
     2322    #skip header
     2323    point_reader.next()
     2324    #read point info from file
     2325    for i,row in enumerate(point_reader):
     2326#        print 'i',i,'row',row
     2327        points.append([float(row[1]),float(row[2])])
     2328        point_name.append(row[0])
     2329       
     2330    #convert to array for file_function
     2331    points_array = array(points,Float)
     2332       
     2333    points_array = ensure_absolute(points_array)
     2334
     2335    dir_name, base = os.path.split(sww_file)   
     2336#    print 'dirname',dir_name, base
     2337    if access(sww_file,R_OK):
     2338        if verbose: print 'File %s exists' %(sww_file)
     2339    else:
     2340        msg = 'File "%s" could not be opened: Error="%s"'\
     2341               %(sww_file, e)
     2342        raise msg
     2343
     2344    sww_files = get_all_swwfiles(look_in_dir=dir_name,
     2345                                 base_name=base,
     2346                                 verbose=verbose)
     2347   
     2348    #to make all the quantities lower case for file_function
     2349    quantities = [quantity.lower() for quantity in quantities]
     2350
     2351    # what is quantities are needed from sww file to calculate output quantities
     2352    # also
     2353
     2354    core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
     2355   
     2356    for sww_file in sww_files:
     2357   
     2358#        print 'sww_file',sww_file
     2359        sww_file = join(dir_name, sww_file+'.sww')
     2360#        print 'sww_file',sww_file, core_quantities
     2361       
     2362        callable_sww = file_function(sww_file,
     2363                                 quantities=core_quantities,
     2364                                 interpolation_points=points_array,
     2365                                 verbose=verbose,
     2366                                 use_cache=use_cache)
     2367
     2368    gauge_file='gauge_'
     2369
     2370    heading = [quantity for quantity in quantities]
     2371    heading.insert(0,'time')
     2372#    print heading, quantities
     2373
     2374    #create a list of csv writers for all the points and write header
     2375    points_writer = []
     2376    for i,point in enumerate(points):
     2377        points_writer.append(writer(file('gauge_'+point_name[i]+'.csv', "wb")))
     2378        points_writer[i].writerow(heading)
     2379
     2380    for time in callable_sww.get_time():
     2381       # points_list = []
     2382
     2383        for point_i, point in enumerate(points_array):
     2384            points_list = [time]
     2385#            print'time',time,'point_i',point_i,point, points_array
     2386            point_quantities = callable_sww(time,point_i)
     2387#            print "quantities", point_quantities
     2388           
     2389#            for i, quantity in enumerate(quantities):
     2390#                points_list.append(quantities)
     2391            for quantity in quantities:
     2392                if quantity==NAN:
     2393                    if verbose: print 'quantity does not exist in' %callable_sww.get_name
     2394                else:
     2395                    if quantity == 'stage':
     2396                        points_list.append(point_quantities[0])
     2397                       
     2398                    if quantity == 'elevation':
     2399                        points_list.append(point_quantities[1])
     2400                       
     2401                    if quantity == 'xmomentum':
     2402                        points_list.append(point_quantities[2])
     2403                       
     2404                    if quantity == 'ymomentum':
     2405                        points_list.append(point_quantities[3])
     2406                       
     2407                    if quantity == 'depth':
     2408                        points_list.append(point_quantities[0] - point_quantities[1])
     2409                       
     2410                       
     2411                    if quantity == 'momentum':
     2412                        momentum = sqrt(point_quantities[2]**2 +\
     2413                                        point_quantities[3]**2)
     2414                        points_list.append(momentum)
     2415                       
     2416                    if quantity == 'speed':
     2417                        #if depth is less than 0.001 then speed = 0.0
     2418                        if point_quantities[0] - point_quantities[1] < 0.001:
     2419                            vel = 0.0
     2420                        else:
     2421                            momentum = sqrt(point_quantities[2]**2 +\
     2422                                            point_quantities[3]**2)
     2423                            vel = momentum/depth             
     2424#                            vel = momentum/(depth + 1.e-6/depth)             
     2425                       
     2426                        points_list.append(vel)
     2427                       
     2428                    if quantity == 'bearing':
     2429                        points_list.append(calc_bearing(point_quantities[2],
     2430                                                        point_quantities[3]))
     2431                       
     2432                #print 'list',points_list
     2433               
     2434            points_writer[point_i].writerow(points_list)
     2435       
     2436
     2437
     2438           
     2439
     2440   
     2441       
Note: See TracChangeset for help on using the changeset viewer.