Changeset 7672


Ignore:
Timestamp:
Mar 26, 2010, 9:26:33 PM (14 years ago)
Author:
hudson
Message:

Wrote 1 failing test for the future gauge_sww2csv centroid intersection option.

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

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

    r7562 r7672  
    15621562       
    15631563
    1564     def test_sww2csv_gauges(self):
    1565 
    1566         def elevation_function(x, y):
    1567             return -x
    1568        
    1569         """Most of this test was copied from test_interpolate
    1570         test_interpole_sww2csv
    1571        
    1572         This is testing the gauge_sww2csv function, by creating a sww file and
    1573         then exporting the gauges and checking the results.
    1574         """
    1575        
    1576         # Create mesh
    1577         mesh_file = tempfile.mktemp(".tsh")   
    1578         points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
    1579         m = Mesh()
    1580         m.add_vertices(points)
    1581         m.auto_segment()
    1582         m.generate_mesh(verbose=False)
    1583         m.export_mesh_file(mesh_file)
    1584        
    1585         # Create shallow water domain
    1586         domain = Domain(mesh_file)
    1587         os.remove(mesh_file)
    1588        
    1589         domain.default_order=2
    1590        
    1591         # This test was made before tight_slope_limiters were introduced
    1592         # Since were are testing interpolation values this is OK
    1593         domain.tight_slope_limiters = 0
    1594        
    1595 
    1596         # Set some field values
    1597         domain.set_quantity('elevation', elevation_function)
    1598         domain.set_quantity('friction', 0.03)
    1599         domain.set_quantity('xmomentum', 3.0)
    1600         domain.set_quantity('ymomentum', 4.0)
    1601 
    1602         ######################
    1603         # Boundary conditions
    1604         B = Transmissive_boundary(domain)
    1605         domain.set_boundary( {'exterior': B})
    1606 
    1607         # This call mangles the stage values.
    1608         domain.distribute_to_vertices_and_edges()
    1609         domain.set_quantity('stage', 1.0)
    1610 
    1611 
    1612         domain.set_name('datatest' + str(time.time()))
    1613         domain.smooth = True
    1614         domain.reduction = mean
    1615 
    1616 
    1617         sww = SWW_file(domain)
    1618         sww.store_connectivity()
    1619         sww.store_timestep()
    1620         domain.set_quantity('stage', 10.0) # This is automatically limited
    1621         # so it will not be less than the elevation
    1622         domain.time = 2.
    1623         sww.store_timestep()
    1624 
    1625         # test the function
    1626         points = [[5.0,1.],[0.5,2.]]
    1627 
    1628         points_file = tempfile.mktemp(".csv")
    1629 #        points_file = 'test_point.csv'
    1630         file_id = open(points_file,"w")
    1631         file_id.write("name, easting, northing, elevation \n\
    1632 point1, 5.0, 1.0, 3.0\n\
    1633 point2, 0.5, 2.0, 9.0\n")
    1634         file_id.close()
    1635 
    1636        
    1637         sww2csv_gauges(sww.filename,
    1638                        points_file,
    1639                        verbose=False,
    1640                        use_cache=False)
    1641 
    1642 #        point1_answers_array = [[0.0,1.0,-5.0,3.0,4.0], [2.0,10.0,-5.0,3.0,4.0]]
    1643         point1_answers_array = [[0.0,0.0,1.0,6.0,-5.0,3.0,4.0], [2.0,2.0/3600.,10.0,15.0,-5.0,3.0,4.0]]
    1644         point1_filename = 'gauge_point1.csv'
    1645         point1_handle = file(point1_filename)
    1646         point1_reader = reader(point1_handle)
    1647         point1_reader.next()
    1648 
    1649         line=[]
    1650         for i,row in enumerate(point1_reader):
    1651 #            print 'i',i,'row',row
    1652             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
    1653                          float(row[4]),float(row[5]),float(row[6])])
    1654 #            print 'assert line',line[i],'point1',point1_answers_array[i]
    1655             assert num.allclose(line[i], point1_answers_array[i])
    1656 
    1657         point2_answers_array = [[0.0,0.0,1.0,1.5,-0.5,3.0,4.0], [2.0,2.0/3600.,10.0,10.5,-0.5,3.0,4.0]]
    1658         point2_filename = 'gauge_point2.csv'
    1659         point2_handle = file(point2_filename)
    1660         point2_reader = reader(point2_handle)
    1661         point2_reader.next()
    1662                        
    1663         line=[]
    1664         for i,row in enumerate(point2_reader):
    1665 #            print 'i',i,'row',row
    1666             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
    1667                          float(row[4]),float(row[5]),float(row[6])])
    1668 #            print 'assert line',line[i],'point1',point1_answers_array[i]
    1669             assert num.allclose(line[i], point2_answers_array[i])
    1670                          
    1671         # clean up
    1672         point1_handle.close()
    1673         point2_handle.close()
    1674         #print "sww.filename",sww.filename
    1675         os.remove(sww.filename)
    1676         os.remove(points_file)
    1677         os.remove(point1_filename)
    1678         os.remove(point2_filename)
    1679 
    1680 
    1681 
    1682     def test_sww2csv_gauges1(self):
    1683         from anuga.pmesh.mesh import Mesh
    1684         from anuga.shallow_water import Domain, Transmissive_boundary
    1685         from csv import reader,writer
    1686         import time
    1687         import string
    1688 
    1689         def elevation_function(x, y):
    1690             return -x
    1691        
    1692         """Most of this test was copied from test_interpolate
    1693         test_interpole_sww2csv
    1694        
    1695         This is testing the gauge_sww2csv function, by creating a sww file and
    1696         then exporting the gauges and checking the results.
    1697        
    1698         This tests the ablity not to have elevation in the points file and
    1699         not store xmomentum and ymomentum
    1700         """
    1701        
    1702         # Create mesh
    1703         mesh_file = tempfile.mktemp(".tsh")   
    1704         points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
    1705         m = Mesh()
    1706         m.add_vertices(points)
    1707         m.auto_segment()
    1708         m.generate_mesh(verbose=False)
    1709         m.export_mesh_file(mesh_file)
    1710        
    1711         # Create shallow water domain
    1712         domain = Domain(mesh_file)
    1713         os.remove(mesh_file)
    1714        
    1715         domain.default_order=2
    1716 
    1717         # Set some field values
    1718         domain.set_quantity('elevation', elevation_function)
    1719         domain.set_quantity('friction', 0.03)
    1720         domain.set_quantity('xmomentum', 3.0)
    1721         domain.set_quantity('ymomentum', 4.0)
    1722 
    1723         ######################
    1724         # Boundary conditions
    1725         B = Transmissive_boundary(domain)
    1726         domain.set_boundary( {'exterior': B})
    1727 
    1728         # This call mangles the stage values.
    1729         domain.distribute_to_vertices_and_edges()
    1730         domain.set_quantity('stage', 1.0)
    1731 
    1732 
    1733         domain.set_name('datatest' + str(time.time()))
    1734         domain.smooth = True
    1735         domain.reduction = mean
    1736 
    1737         sww = SWW_file(domain)
    1738         sww.store_connectivity()
    1739         sww.store_timestep()
    1740         domain.set_quantity('stage', 10.0) # This is automatically limited
    1741         # so it will not be less than the elevation
    1742         domain.time = 2.
    1743         sww.store_timestep()
    1744 
    1745         # test the function
    1746         points = [[5.0,1.],[0.5,2.]]
    1747 
    1748         points_file = tempfile.mktemp(".csv")
    1749 #        points_file = 'test_point.csv'
    1750         file_id = open(points_file,"w")
    1751         file_id.write("name,easting,northing \n\
    1752 point1, 5.0, 1.0\n\
    1753 point2, 0.5, 2.0\n")
    1754         file_id.close()
    1755 
    1756         sww2csv_gauges(sww.filename,
    1757                             points_file,
    1758                             quantities=['stage', 'elevation'],
    1759                             use_cache=False,
    1760                             verbose=False)
    1761 
    1762         point1_answers_array = [[0.0,1.0,-5.0], [2.0,10.0,-5.0]]
    1763         point1_filename = 'gauge_point1.csv'
    1764         point1_handle = file(point1_filename)
    1765         point1_reader = reader(point1_handle)
    1766         point1_reader.next()
    1767 
    1768         line=[]
    1769         for i,row in enumerate(point1_reader):
    1770 #            print 'i',i,'row',row
    1771             # note the 'hole' (element 1) below - skip the new 'hours' field
    1772             line.append([float(row[0]),float(row[2]),float(row[3])])
    1773             #print 'line',line[i],'point1',point1_answers_array[i]
    1774             assert num.allclose(line[i], point1_answers_array[i])
    1775 
    1776         point2_answers_array = [[0.0,1.0,-0.5], [2.0,10.0,-0.5]]
    1777         point2_filename = 'gauge_point2.csv'
    1778         point2_handle = file(point2_filename)
    1779         point2_reader = reader(point2_handle)
    1780         point2_reader.next()
    1781                        
    1782         line=[]
    1783         for i,row in enumerate(point2_reader):
    1784 #            print 'i',i,'row',row
    1785             # note the 'hole' (element 1) below - skip the new 'hours' field
    1786             line.append([float(row[0]),float(row[2]),float(row[3])])
    1787 #            print 'line',line[i],'point1',point1_answers_array[i]
    1788             assert num.allclose(line[i], point2_answers_array[i])
    1789                          
    1790         # clean up
    1791         point1_handle.close()
    1792         point2_handle.close()
    1793         #print "sww.filename",sww.filename
    1794         os.remove(sww.filename)
    1795         os.remove(points_file)
    1796         os.remove(point1_filename)
    1797         os.remove(point2_filename)
    1798 
    1799 
    1800     def test_sww2csv_gauges2(self):
    1801 
    1802         def elevation_function(x, y):
    1803             return -x
    1804        
    1805         """Most of this test was copied from test_interpolate
    1806         test_interpole_sww2csv
    1807        
    1808         This is testing the gauge_sww2csv function, by creating a sww file and
    1809         then exporting the gauges and checking the results.
    1810        
    1811         This is the same as sww2csv_gauges except set domain.set_starttime to 5.
    1812         Therefore testing the storing of the absolute time in the csv files
    1813         """
    1814        
    1815         # Create mesh
    1816         mesh_file = tempfile.mktemp(".tsh")   
    1817         points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
    1818         m = Mesh()
    1819         m.add_vertices(points)
    1820         m.auto_segment()
    1821         m.generate_mesh(verbose=False)
    1822         m.export_mesh_file(mesh_file)
    1823        
    1824         # Create shallow water domain
    1825         domain = Domain(mesh_file)
    1826         os.remove(mesh_file)
    1827        
    1828         domain.default_order=2
    1829 
    1830         # This test was made before tight_slope_limiters were introduced
    1831         # Since were are testing interpolation values this is OK
    1832         domain.tight_slope_limiters = 0         
    1833 
    1834         # Set some field values
    1835         domain.set_quantity('elevation', elevation_function)
    1836         domain.set_quantity('friction', 0.03)
    1837         domain.set_quantity('xmomentum', 3.0)
    1838         domain.set_quantity('ymomentum', 4.0)
    1839         domain.set_starttime(5)
    1840 
    1841         ######################
    1842         # Boundary conditions
    1843         B = Transmissive_boundary(domain)
    1844         domain.set_boundary( {'exterior': B})
    1845 
    1846         # This call mangles the stage values.
    1847         domain.distribute_to_vertices_and_edges()
    1848         domain.set_quantity('stage', 1.0)
    1849        
    1850 
    1851 
    1852         domain.set_name('datatest' + str(time.time()))
    1853         domain.smooth = True
    1854         domain.reduction = mean
    1855 
    1856         sww = SWW_file(domain)
    1857         sww.store_connectivity()
    1858         sww.store_timestep()
    1859         domain.set_quantity('stage', 10.0) # This is automatically limited
    1860         # so it will not be less than the elevation
    1861         domain.time = 2.
    1862         sww.store_timestep()
    1863 
    1864         # test the function
    1865         points = [[5.0,1.],[0.5,2.]]
    1866 
    1867         points_file = tempfile.mktemp(".csv")
    1868 #        points_file = 'test_point.csv'
    1869         file_id = open(points_file,"w")
    1870         file_id.write("name, easting, northing, elevation \n\
    1871 point1, 5.0, 1.0, 3.0\n\
    1872 point2, 0.5, 2.0, 9.0\n")
    1873         file_id.close()
    1874 
    1875        
    1876         sww2csv_gauges(sww.filename,
    1877                             points_file,
    1878                             verbose=False,
    1879                             use_cache=False)
    1880 
    1881 #        point1_answers_array = [[0.0,1.0,-5.0,3.0,4.0], [2.0,10.0,-5.0,3.0,4.0]]
    1882         point1_answers_array = [[5.0,5.0/3600.,1.0,6.0,-5.0,3.0,4.0], [7.0,7.0/3600.,10.0,15.0,-5.0,3.0,4.0]]
    1883         point1_filename = 'gauge_point1.csv'
    1884         point1_handle = file(point1_filename)
    1885         point1_reader = reader(point1_handle)
    1886         point1_reader.next()
    1887 
    1888         line=[]
    1889         for i,row in enumerate(point1_reader):
    1890             #print 'i',i,'row',row
    1891             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
    1892                          float(row[4]), float(row[5]), float(row[6])])
    1893             #print 'assert line',line[i],'point1',point1_answers_array[i]
    1894             assert num.allclose(line[i], point1_answers_array[i])
    1895 
    1896         point2_answers_array = [[5.0,5.0/3600.,1.0,1.5,-0.5,3.0,4.0], [7.0,7.0/3600.,10.0,10.5,-0.5,3.0,4.0]]
    1897         point2_filename = 'gauge_point2.csv'
    1898         point2_handle = file(point2_filename)
    1899         point2_reader = reader(point2_handle)
    1900         point2_reader.next()
    1901                        
    1902         line=[]
    1903         for i,row in enumerate(point2_reader):
    1904             #print 'i',i,'row',row
    1905             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
    1906                          float(row[4]),float(row[5]), float(row[6])])
    1907             #print 'assert line',line[i],'point1',point1_answers_array[i]
    1908             assert num.allclose(line[i], point2_answers_array[i])
    1909                          
    1910         # clean up
    1911         point1_handle.close()
    1912         point2_handle.close()
    1913         #print "sww.filename",sww.filename
    1914         os.remove(sww.filename)
    1915         os.remove(points_file)
    1916         os.remove(point1_filename)
    1917         os.remove(point2_filename)
    1918 
    19191564
    19201565    def test_greens_law(self):

  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r7317 r7672  
    4343# @param use_cache True means that caching of intermediate result is attempted.
    4444# @param boundary_polygon
     45# @param output_centroids if True, data for the centroid of the triangle will be output
    4546# @return A callable object.
    4647def file_function(filename,
    … …  
    5253                  verbose=False,
    5354                  use_cache=False,
    54                   boundary_polygon=None):
     55                  boundary_polygon=None,
     56                                  output_centroids=False):
    5557    """Read time history of spatial data from NetCDF file and return
    5658    a callable object.
    … …  
    137139              'time_limit': time_limit,                                 
    138140              'verbose': verbose,
    139               'boundary_polygon': boundary_polygon}
     141              'boundary_polygon': boundary_polygon,
     142                          'output_centroids': output_centroids}
    140143
    141144    # Call underlying engine with or without caching
    … …  
    200203                   time_limit=None,
    201204                   verbose=False,
    202                    boundary_polygon=None):
     205                   boundary_polygon=None,
     206                                   output_centroids=False):
    203207    """Internal function
    204208   
    … …  
    227231                                        time_limit=time_limit,
    228232                                        verbose=verbose,
    229                                         boundary_polygon=boundary_polygon)
     233                                        boundary_polygon=boundary_polygon,
     234                                                                                output_centroids=output_centroids)
    230235    else:
    231236        # FIXME (Ole): Could add csv file here to address Ted Rigby's
    … …  
    254259                             time_limit=None,           
    255260                             verbose=False,
    256                              boundary_polygon=None):
     261                             boundary_polygon=None,
     262                                                         output_centroids=False):
    257263    """Read time history of spatial data from NetCDF sww file and
    258264    return a callable object f(t,x,y)
    … …  
    739745##
    740746# @brief Read a .sww file and plot the time series.
    741 # @param swwfiles Dictionary of .sww files.
    742 # @param gauge_filename Name of gauge file.
    743 # @param production_dirs ??
    744 # @param report If True, write figures to report directory.
    745 # @param reportname Name of generated report (if any).
    746 # @param plot_quantity List containing quantities to plot.
    747 # @param generate_fig If True, generate figures as well as CSV files.
    748 # @param surface If True, then generate solution surface with 3d plot.
    749 # @param time_min Beginning of user defined time range for plotting purposes.
    750 # @param time_max End of user defined time range for plotting purposes.
    751 # @param time_thinning ??
    752 # @param time_unit ??
    753 # @param title_on If True, export standard graphics with title.
    754 # @param use_cache If True, use caching.
    755 # @param verbose If True, this function is verbose.
     747# @note This function is deprecated - use gauge.sww2timeseries instead.
     748#
    756749def sww2timeseries(swwfiles,
    757750                   gauge_filename,
    … …  
    769762                   use_cache=False,
    770763                   verbose=False):
    771     """ Read sww file and plot the time series for the
    772     prescribed quantities at defined gauge locations and
    773     prescribed time range.
    774 
    775     Input variables:
    776 
    777     swwfiles        - dictionary of sww files with label_ids (used in
    778                       generating latex output. It will be part of
    779                       the directory name of file_loc (typically the timestamp).
    780                       Helps to differentiate latex files for different
    781                       simulations for a particular scenario. 
    782                     - assume that all conserved quantities have been stored
    783                     - assume each sww file has been simulated with same timestep
    784    
    785     gauge_filename  - name of file containing gauge data
    786                         - easting, northing, name , elevation?
    787                     - OR (this is not yet done)
    788                         - structure which can be converted to a numeric array,
    789                           such as a geospatial data object
    790                      
    791     production_dirs -  A list of list, example {20061101_121212: '1 in 10000',
    792                                                 'boundaries': 'urs boundary'}
    793                       this will use the second part as the label and the
    794                       first part as the ?
    795                       #FIXME: Is it a list or a dictionary
    796                       # This is probably obsolete by now
    797                      
    798     report          - if True, then write figures to report_figures directory in
    799                       relevant production directory
    800                     - if False, figures are already stored with sww file
    801                     - default to False
    802 
    803     reportname      - name for report if wishing to generate report
    804    
    805     plot_quantity   - list containing quantities to plot, they must
    806                       be the name of an existing quantity or one of
    807                       the following possibilities
    808                     - possibilities:
    809                         - stage; 'stage'
    810                         - depth; 'depth'
    811                         - speed; calculated as absolute momentum
    812                          (pointwise) divided by depth; 'speed'
    813                         - bearing; calculated as the angle of the momentum
    814                           vector (xmomentum, ymomentum) from the North; 'bearing'
    815                         - absolute momentum; calculated as
    816                           sqrt(xmomentum^2 + ymomentum^2); 'momentum'
    817                         - x momentum; 'xmomentum'
    818                         - y momentum; 'ymomentum'
    819                     - default will be ['stage', 'speed', 'bearing']
    820 
    821     generate_fig     - if True, generate figures as well as csv file
    822                      - if False, csv files created only
    823                      
    824     surface          - if True, then generate solution surface with 3d plot
    825                        and save to current working directory
    826                      - default = False
    827    
    828     time_min         - beginning of user defined time range for plotting purposes
    829                         - default will be first available time found in swwfile
    830                        
    831     time_max         - end of user defined time range for plotting purposes
    832                         - default will be last available time found in swwfile
    833                        
    834     title_on        - if True, export standard graphics with title
    835                     - if False, export standard graphics without title
    836 
    837 
    838     Output:
    839    
    840     - time series data stored in .csv for later use if required.
    841       Name = gauges_timeseries followed by gauge name
    842     - latex file will be generated in same directory as where script is
    843       run (usually production scenario directory.
    844       Name = latexoutputlabel_id.tex
    845 
    846     Other important information:
    847    
    848     It is assumed that the used has stored all the conserved quantities
    849     and elevation during the scenario run, i.e.
    850     ['stage', 'elevation', 'xmomentum', 'ymomentum']
    851     If this has not occurred then sww2timeseries will not work.
    852 
    853 
    854     Usage example
    855     texname = sww2timeseries({project.boundary_name + '.sww': ''},
    856                              project.polygons_dir + sep + 'boundary_extent.csv',
    857                              project.anuga_dir,
    858                              report = False,
    859                              plot_quantity = ['stage', 'speed', 'bearing'],
    860                              time_min = None,
    861                              time_max = None,
    862                              title_on = True,   
    863                              verbose = True)
    864    
    865     """
    866 
    867     msg = 'NOTE: A new function is available to create csv files from sww '
    868     msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util'
    869     msg += ' PLUS another new function to create graphs from csv files called '
    870     msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util'
    871     log.critical(msg)
    872    
    873     k = _sww2timeseries(swwfiles,
    874                         gauge_filename,
    875                         production_dirs,
    876                         report,
    877                         reportname,
    878                         plot_quantity,
    879                         generate_fig,
    880                         surface,
    881                         time_min,
    882                         time_max,
    883                         time_thinning,                       
    884                         time_unit,
    885                         title_on,
    886                         use_cache,
    887                         verbose)
    888     return k
    889 
    890 
    891 ##
    892 # @brief Read a .sww file and plot the time series.
    893 # @param swwfiles Dictionary of .sww files.
    894 # @param gauge_filename Name of gauge file.
    895 # @param production_dirs ??
    896 # @param report If True, write figures to report directory.
    897 # @param reportname Name of generated report (if any).
    898 # @param plot_quantity List containing quantities to plot.
    899 # @param generate_fig If True, generate figures as well as CSV files.
    900 # @param surface If True, then generate solution surface with 3d plot.
    901 # @param time_min Beginning of user defined time range for plotting purposes.
    902 # @param time_max End of user defined time range for plotting purposes.
    903 # @param time_thinning ??
    904 # @param time_unit ??
    905 # @param title_on If True, export standard graphics with title.
    906 # @param use_cache If True, use caching.
    907 # @param verbose If True, this function is verbose.
    908 def _sww2timeseries(swwfiles,
    909                     gauge_filename,
    910                     production_dirs,
    911                     report = None,
    912                     reportname = None,
    913                     plot_quantity = None,
    914                     generate_fig = False,
    915                     surface = None,
    916                     time_min = None,
    917                     time_max = None,
    918                     time_thinning = 1,                   
    919                     time_unit = None,
    920                     title_on = None,
    921                     use_cache = False,
    922                     verbose = False):   
    923        
    924     # FIXME(Ole): Shouldn't print statements here be governed by verbose?
    925     assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
    926    
    927     try:
    928         fid = open(gauge_filename)
    929     except Exception, e:
    930         msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e)
    931         raise msg
    932 
    933     if report is None:
    934         report = False
    935        
    936     if plot_quantity is None:
    937         plot_quantity = ['depth', 'speed']
    938     else:
    939         assert type(plot_quantity) == list, 'plot_quantity must be a list'
    940         check_list(plot_quantity)
    941 
    942     if surface is None:
    943         surface = False
    944 
    945     if time_unit is None:
    946         time_unit = 'hours'
    947    
    948     if title_on is None:
    949         title_on = True
    950    
    951     if verbose: log.critical('Gauges obtained from: %s' % gauge_filename)
    952 
    953     gauges, locations, elev = get_gauges_from_file(gauge_filename)
    954 
    955     sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
    956 
    957     file_loc = []
    958     f_list = []
    959     label_id = []
    960     leg_label = []
    961     themaxT = 0.0
    962     theminT = 0.0
    963 
    964     for swwfile in swwfiles.keys():
    965         try:
    966             fid = open(swwfile)
    967         except Exception, e:
    968             msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e)
    969             raise msg
    970 
    971         if verbose:
    972             log.critical('swwfile = %s' % swwfile)
    973 
    974         # Extract parent dir name and use as label
    975         path, _ = os.path.split(swwfile)
    976         _, label = os.path.split(path)       
    977        
    978         leg_label.append(label)
    979 
    980         f = file_function(swwfile,
    981                           quantities = sww_quantity,
    982                           interpolation_points = gauges,
    983                           time_thinning = time_thinning,
    984                           verbose = verbose,
    985                           use_cache = use_cache)
    986 
    987         # determine which gauges are contained in sww file
    988         count = 0
    989         gauge_index = []
    990         for k, g in enumerate(gauges):
    991             if f(0.0, point_id = k)[2] > 1.0e6:
    992                 count += 1
    993                 if count == 1: log.critical('Gauges not contained here:')
    994                 log.critical(locations[k])
    995             else:
    996                 gauge_index.append(k)
    997 
    998         if len(gauge_index) > 0:
    999             log.critical('Gauges contained here:')
    1000         else:
    1001             log.critical('No gauges contained here.')
    1002         for i in range(len(gauge_index)):
    1003              log.critical(locations[gauge_index[i]])
    1004              
    1005         index = swwfile.rfind(sep)
    1006         file_loc.append(swwfile[:index+1])
    1007         label_id.append(swwfiles[swwfile])
    1008        
    1009         f_list.append(f)
    1010         maxT = max(f.get_time())
    1011         minT = min(f.get_time())
    1012         if maxT > themaxT: themaxT = maxT
    1013         if minT > theminT: theminT = minT
    1014 
    1015     if time_min is None:
    1016         time_min = theminT # min(T)
    1017     else:
    1018         if time_min < theminT: # min(T):
    1019             msg = 'Minimum time entered not correct - please try again'
    1020             raise Exception, msg
    1021 
    1022     if time_max is None:
    1023         time_max = themaxT # max(T)
    1024     else:
    1025         if time_max > themaxT: # max(T):
    1026             msg = 'Maximum time entered not correct - please try again'
    1027             raise Exception, msg
    1028 
    1029     if verbose and len(gauge_index) > 0:
    1030          log.critical('Inputs OK - going to generate figures')
    1031 
    1032     if len(gauge_index) <> 0:
    1033         texfile, elev_output = \
    1034             generate_figures(plot_quantity, file_loc, report, reportname,
    1035                              surface, leg_label, f_list, gauges, locations,
    1036                              elev, gauge_index, production_dirs, time_min,
    1037                              time_max, time_unit, title_on, label_id,
    1038                              generate_fig, verbose)
    1039     else:
    1040         texfile = ''
    1041         elev_output = []
    1042 
    1043     return texfile, elev_output
     764        return gauge.sww2timeseries(swwfiles, gauge_filename, production_dirs, report, reportname, \
     765                                    plot_quantity, generate_fig, surface, time_min, time_max,     \
     766                                                                time_thinning, time_unit, title_on, use_cache, verbose)
     767
    1044768
    1045769
    … …  
    1049773# @return A (gauges, gaugelocation, elev) tuple.
    1050774def get_gauges_from_file(filename):
    1051     """ Read in gauge information from file
    1052     """
    1053 
    1054     from os import sep, getcwd, access, F_OK, mkdir
    1055 
    1056     # Get data from the gauge file
    1057     fid = open(filename)
    1058     lines = fid.readlines()
    1059     fid.close()
    1060    
    1061     gauges = []
    1062     gaugelocation = []
    1063     elev = []
    1064 
    1065     # Check header information   
    1066     line1 = lines[0]
    1067     line11 = line1.split(',')
    1068 
    1069     if isinstance(line11[0], str) is True:
    1070         # We have found text in the first line
    1071         east_index = None
    1072         north_index = None
    1073         name_index = None
    1074         elev_index = None
    1075 
    1076         for i in range(len(line11)):
    1077             if line11[i].strip().lower() == 'easting':   east_index = i
    1078             if line11[i].strip().lower() == 'northing':  north_index = i
    1079             if line11[i].strip().lower() == 'name':      name_index = i
    1080             if line11[i].strip().lower() == 'elevation': elev_index = i
    1081 
    1082         if east_index < len(line11) and north_index < len(line11):
    1083             pass
    1084         else:
    1085             msg = 'WARNING: %s does not contain correct header information' \
    1086                   % filename
    1087             msg += 'The header must be: easting, northing, name, elevation'
    1088             raise Exception, msg
    1089 
    1090         if elev_index is None:
    1091             raise Exception
    1092    
    1093         if name_index is None:
    1094             raise Exception
    1095 
    1096         lines = lines[1:] # Remove header from data
    1097     else:
    1098         # No header, assume that this is a simple easting, northing file
    1099 
    1100         msg = 'There was no header in file %s and the number of columns is %d' \
    1101               % (filename, len(line11))
    1102         msg += '- was assuming two columns corresponding to Easting and Northing'
    1103         assert len(line11) == 2, msg
    1104 
    1105         east_index = 0
    1106         north_index = 1
    1107 
    1108         N = len(lines)
    1109         elev = [-9999]*N
    1110         gaugelocation = range(N)
    1111        
    1112     # Read in gauge data
    1113     for line in lines:
    1114         fields = line.split(',')
    1115 
    1116         gauges.append([float(fields[east_index]), float(fields[north_index])])
    1117 
    1118         if len(fields) > 2:
    1119             elev.append(float(fields[elev_index]))
    1120             loc = fields[name_index]
    1121             gaugelocation.append(loc.strip('\n'))
    1122 
    1123     return gauges, gaugelocation, elev
     775    return gauge.get_from_file(filename)
    1124776
    1125777
    … …  
    24252077    This is really returning speed, not velocity.
    24262078    """
    2427    
    2428     from csv import reader,writer
    2429     from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN
    2430     import string
    2431     from anuga.shallow_water.data_manager import get_all_swwfiles
    2432 
    2433     assert type(gauge_file) == type(''), 'Gauge filename must be a string'
    2434     assert type(out_name) == type(''), 'Output filename prefix must be a string'
    2435    
    2436     try:
    2437         point_reader = reader(file(gauge_file))
    2438     except Exception, e:
    2439         msg = 'File "%s" could not be opened: Error="%s"' % (gauge_file, e)
    2440         raise msg
    2441 
    2442     if verbose: log.critical('Gauges obtained from: %s' % gauge_file)
    2443    
    2444     point_reader = reader(file(gauge_file))
    2445     points = []
    2446     point_name = []
    2447    
    2448     # read point info from file
    2449     for i,row in enumerate(point_reader):
    2450         # read header and determine the column numbers to read correctly.
    2451         if i==0:
    2452             for j,value in enumerate(row):
    2453                 if value.strip()=='easting':easting=j
    2454                 if value.strip()=='northing':northing=j
    2455                 if value.strip()=='name':name=j
    2456                 if value.strip()=='elevation':elevation=j
    2457         else:
    2458             points.append([float(row[easting]),float(row[northing])])
    2459             point_name.append(row[name])
    2460        
    2461     #convert to array for file_function
    2462     points_array = num.array(points,num.float)
    2463        
    2464     points_array = ensure_absolute(points_array)
    2465 
    2466     dir_name, base = os.path.split(sww_file)   
    2467 
    2468     #need to get current directory so when path and file
    2469     #are "joined" below the directory is correct
    2470     if dir_name == '':
    2471         dir_name =getcwd()
    2472        
    2473     if access(sww_file,R_OK):
    2474         if verbose: log.critical('File %s exists' % sww_file)
    2475     else:
    2476         msg = 'File "%s" could not be opened: no read permission' % sww_file
    2477         raise msg
    2478 
    2479     sww_files = get_all_swwfiles(look_in_dir=dir_name,
    2480                                  base_name=base,
    2481                                  verbose=verbose)
    2482 
    2483     # fudge to get SWW files in 'correct' order, oldest on the left
    2484     sww_files.sort()
    2485 
    2486     if verbose:
    2487         log.critical('sww files=%s' % sww_files)
    2488    
    2489     #to make all the quantities lower case for file_function
    2490     quantities = [quantity.lower() for quantity in quantities]
    2491 
    2492     # what is quantities are needed from sww file to calculate output quantities
    2493     # also
    2494 
    2495     core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
    2496 
    2497     gauge_file = out_name
    2498 
    2499     heading = [quantity for quantity in quantities]
    2500     heading.insert(0,'time')
    2501     heading.insert(1,'hours')
    2502 
    2503     #create a list of csv writers for all the points and write header
    2504     points_writer = []
    2505     for point_i,point in enumerate(points):
    2506         points_writer.append(writer(file(dir_name + sep + gauge_file
    2507                                          + point_name[point_i] + '.csv', "wb")))
    2508         points_writer[point_i].writerow(heading)
    2509    
    2510     if verbose: log.critical('Writing csv files')
    2511 
    2512     quake_offset_time = None
    2513 
    2514     for sww_file in sww_files:
    2515         sww_file = join(dir_name, sww_file+'.sww')
    2516         callable_sww = file_function(sww_file,
    2517                                      quantities=core_quantities,
    2518                                      interpolation_points=points_array,
    2519                                      verbose=verbose,
    2520                                      use_cache=use_cache)
    2521 
    2522         if quake_offset_time is None:
    2523             quake_offset_time = callable_sww.starttime
    2524 
    2525         for time in callable_sww.get_time():
    2526             for point_i, point in enumerate(points_array):
    2527                 #add domain starttime to relative time.
    2528                 quake_time = time + quake_offset_time
    2529                 points_list = [quake_time, quake_time/3600.]# fudge around SWW time bug
    2530                 point_quantities = callable_sww(time,point_i)
    2531                
    2532                 for quantity in quantities:
    2533                     if quantity == NAN:
    2534                         log.critical('quantity does not exist in %s'
    2535                                      % callable_sww.get_name)
    2536                     else:
    2537                         if quantity == 'stage':
    2538                             points_list.append(point_quantities[0])
    2539                            
    2540                         if quantity == 'elevation':
    2541                             points_list.append(point_quantities[1])
    2542                            
    2543                         if quantity == 'xmomentum':
    2544                             points_list.append(point_quantities[2])
    2545                            
    2546                         if quantity == 'ymomentum':
    2547                             points_list.append(point_quantities[3])
    2548                            
    2549                         if quantity == 'depth':
    2550                             points_list.append(point_quantities[0]
    2551                                                - point_quantities[1])
    2552 
    2553                         if quantity == 'momentum':
    2554                             momentum = sqrt(point_quantities[2]**2
    2555                                             + point_quantities[3]**2)
    2556                             points_list.append(momentum)
    2557                            
    2558                         if quantity == 'speed':
    2559                             #if depth is less than 0.001 then speed = 0.0
    2560                             if point_quantities[0] - point_quantities[1] < 0.001:
    2561                                 vel = 0.0
    2562                             else:
    2563                                 if point_quantities[2] < 1.0e6:
    2564                                     momentum = sqrt(point_quantities[2]**2
    2565                                                     + point_quantities[3]**2)
    2566                                     vel = momentum / (point_quantities[0]
    2567                                                       - point_quantities[1])
    2568                                 else:
    2569                                     momentum = 0
    2570                                     vel = 0
    2571                                
    2572                             points_list.append(vel)
    2573                            
    2574                         if quantity == 'bearing':
    2575                             points_list.append(calc_bearing(point_quantities[2],
    2576                                                             point_quantities[3]))
    2577 
    2578                 points_writer[point_i].writerow(points_list)
    2579 
     2079    from gauge import sww2csv
     2080       
     2081    sww2csv(sww_file, gauge_file, out_name, quantities, verbose, use_cache)
     2082       
     2083       
    25802084##
    25812085# @brief Get a wave height at a certain depth given wave height at another depth.
Note: See TracChangeset for help on using the changeset viewer.