Changeset 7672
- Timestamp:
- Mar 26, 2010, 9:26:33 PM (15 years ago)
- Location:
- anuga_core/source/anuga/abstract_2d_finite_volumes
- Files:
-
- 2 added
- 2 edited
-
gauge.py (added)
-
test_gauge.py (added)
-
test_util.py (modified) (1 diff)
-
util.py (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py
r7562 r7672 1562 1562 1563 1563 1564 def test_sww2csv_gauges(self):1565 1566 def elevation_function(x, y):1567 return -x1568 1569 """Most of this test was copied from test_interpolate1570 test_interpole_sww2csv1571 1572 This is testing the gauge_sww2csv function, by creating a sww file and1573 then exporting the gauges and checking the results.1574 """1575 1576 # Create mesh1577 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 domain1586 domain = Domain(mesh_file)1587 os.remove(mesh_file)1588 1589 domain.default_order=21590 1591 # This test was made before tight_slope_limiters were introduced1592 # Since were are testing interpolation values this is OK1593 domain.tight_slope_limiters = 01594 1595 1596 # Set some field values1597 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 conditions1604 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 = True1614 domain.reduction = mean1615 1616 1617 sww = SWW_file(domain)1618 sww.store_connectivity()1619 sww.store_timestep()1620 domain.set_quantity('stage', 10.0) # This is automatically limited1621 # so it will not be less than the elevation1622 domain.time = 2.1623 sww.store_timestep()1624 1625 # test the function1626 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',row1652 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',row1666 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 up1672 point1_handle.close()1673 point2_handle.close()1674 #print "sww.filename",sww.filename1675 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 Mesh1684 from anuga.shallow_water import Domain, Transmissive_boundary1685 from csv import reader,writer1686 import time1687 import string1688 1689 def elevation_function(x, y):1690 return -x1691 1692 """Most of this test was copied from test_interpolate1693 test_interpole_sww2csv1694 1695 This is testing the gauge_sww2csv function, by creating a sww file and1696 then exporting the gauges and checking the results.1697 1698 This tests the ablity not to have elevation in the points file and1699 not store xmomentum and ymomentum1700 """1701 1702 # Create mesh1703 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 domain1712 domain = Domain(mesh_file)1713 os.remove(mesh_file)1714 1715 domain.default_order=21716 1717 # Set some field values1718 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 conditions1725 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 = True1735 domain.reduction = mean1736 1737 sww = SWW_file(domain)1738 sww.store_connectivity()1739 sww.store_timestep()1740 domain.set_quantity('stage', 10.0) # This is automatically limited1741 # so it will not be less than the elevation1742 domain.time = 2.1743 sww.store_timestep()1744 1745 # test the function1746 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',row1771 # note the 'hole' (element 1) below - skip the new 'hours' field1772 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',row1785 # note the 'hole' (element 1) below - skip the new 'hours' field1786 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 up1791 point1_handle.close()1792 point2_handle.close()1793 #print "sww.filename",sww.filename1794 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 -x1804 1805 """Most of this test was copied from test_interpolate1806 test_interpole_sww2csv1807 1808 This is testing the gauge_sww2csv function, by creating a sww file and1809 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 files1813 """1814 1815 # Create mesh1816 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 domain1825 domain = Domain(mesh_file)1826 os.remove(mesh_file)1827 1828 domain.default_order=21829 1830 # This test was made before tight_slope_limiters were introduced1831 # Since were are testing interpolation values this is OK1832 domain.tight_slope_limiters = 01833 1834 # Set some field values1835 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 conditions1843 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 = True1854 domain.reduction = mean1855 1856 sww = SWW_file(domain)1857 sww.store_connectivity()1858 sww.store_timestep()1859 domain.set_quantity('stage', 10.0) # This is automatically limited1860 # so it will not be less than the elevation1861 domain.time = 2.1862 sww.store_timestep()1863 1864 # test the function1865 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',row1891 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',row1905 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 up1911 point1_handle.close()1912 point2_handle.close()1913 #print "sww.filename",sww.filename1914 os.remove(sww.filename)1915 os.remove(points_file)1916 os.remove(point1_filename)1917 os.remove(point2_filename)1918 1919 1564 1920 1565 def test_greens_law(self): -
anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
r7317 r7672 43 43 # @param use_cache True means that caching of intermediate result is attempted. 44 44 # @param boundary_polygon 45 # @param output_centroids if True, data for the centroid of the triangle will be output 45 46 # @return A callable object. 46 47 def file_function(filename, … … 52 53 verbose=False, 53 54 use_cache=False, 54 boundary_polygon=None): 55 boundary_polygon=None, 56 output_centroids=False): 55 57 """Read time history of spatial data from NetCDF file and return 56 58 a callable object. … … 137 139 'time_limit': time_limit, 138 140 'verbose': verbose, 139 'boundary_polygon': boundary_polygon} 141 'boundary_polygon': boundary_polygon, 142 'output_centroids': output_centroids} 140 143 141 144 # Call underlying engine with or without caching … … 200 203 time_limit=None, 201 204 verbose=False, 202 boundary_polygon=None): 205 boundary_polygon=None, 206 output_centroids=False): 203 207 """Internal function 204 208 … … 227 231 time_limit=time_limit, 228 232 verbose=verbose, 229 boundary_polygon=boundary_polygon) 233 boundary_polygon=boundary_polygon, 234 output_centroids=output_centroids) 230 235 else: 231 236 # FIXME (Ole): Could add csv file here to address Ted Rigby's … … 254 259 time_limit=None, 255 260 verbose=False, 256 boundary_polygon=None): 261 boundary_polygon=None, 262 output_centroids=False): 257 263 """Read time history of spatial data from NetCDF sww file and 258 264 return a callable object f(t,x,y) … … 739 745 ## 740 746 # @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 # 756 749 def sww2timeseries(swwfiles, 757 750 gauge_filename, … … 769 762 use_cache=False, 770 763 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 1044 768 1045 769 … … 1049 773 # @return A (gauges, gaugelocation, elev) tuple. 1050 774 def 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) 1124 776 1125 777 … … 2425 2077 This is really returning speed, not velocity. 2426 2078 """ 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 2580 2084 ## 2581 2085 # @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.
