Changeset 4876
- Timestamp:
- Dec 6, 2007, 5:37:13 PM (17 years ago)
- 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 15 15 16 16 from sys import platform 17 18 from anuga.pmesh.mesh import Mesh 19 from anuga.shallow_water import Domain, Transmissive_boundary 20 from anuga.shallow_water.data_manager import get_dataobject 21 from csv import reader,writer 22 import time 23 import string 17 24 18 25 def test_function(x, y): … … 1392 1399 1393 1400 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\ 1464 point1, 5.0, 1.0, 3.0\n\ 1465 point2, 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\ 1583 point1, 5.0, 1.0\n\ 1584 point2, 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) 1395 1627 1396 1628 #------------------------------------------------------------- 1397 1629 if __name__ == "__main__": 1398 1630 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) 1401 1634 runner.run(suite) 1402 1635 -
anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
r4820 r4876 9 9 import os 10 10 11 from os import remove, mkdir, access, F_OK, W_OK, sep,mkdir11 from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep,mkdir 12 12 from os.path import exists, basename, split,join 13 13 from warnings import warn … … 1494 1494 return xc 1495 1495 1496 def make_plots_from_csv_file(directories_dic={},1496 def csv2timeseries_graphs(directories_dic={}, 1497 1497 output_dir='', 1498 1498 base_name=None, 1499 1499 plot_numbers='0', 1500 quantities=[' Stage'],1500 quantities=['stage'], 1501 1501 extra_plot_name='', 1502 1502 assess_all_csv_files=True, … … 1504 1504 verbose=False): 1505 1505 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 1507 1511 plot time series such as Stage, Speed, etc. Will also plot several 1508 1512 time series on one plot. Filenames must follow this convention, … … 1565 1569 #be very good and very flexable.... but this works for now! 1566 1570 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... 1588 1573 1589 1574 seconds_in_hour = 3600 … … 1592 1577 speed_label = 'speed (m/s)' 1593 1578 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)' 1594 1583 1595 1584 if extra_plot_name != '': … … 1606 1595 1607 1596 #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. 1609 1599 max_start_time= 0. 1610 1600 min_start_time = 100000 … … 1647 1637 directory_add_tide = directories_dic[directory][2] 1648 1638 1649 time = [float(x) for x in attribute_dic[" Time"]]1639 time = [float(x) for x in attribute_dic["time"]] 1650 1640 min_t, max_t = get_min_max_values(time,min_t,max_t) 1651 1641 1652 stage = [float(x) for x in attribute_dic[" Stage"]]1642 stage = [float(x) for x in attribute_dic["stage"]] 1653 1643 stage =array(stage)+directory_add_tide 1654 1644 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) 1655 1649 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, 1661 1653 min_mom, 1662 1654 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) 1663 1675 1664 1676 # print 'min_sp, max_sp',min_sp, max_sp, … … 1676 1688 (max_t+max_start_time)/seconds_in_hour, 1677 1689 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) 1678 1696 momentum_axis = (min_start_time/seconds_in_hour, 1679 1697 (max_t+max_start_time)/seconds_in_hour, 1680 1698 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 1685 1706 cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k'] 1686 1707 … … 1704 1725 attribute_dic, title_index_dic = csv2dict(file+'.csv') 1705 1726 #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"]] 1707 1728 1708 1729 #do maths to list by changing to array 1709 1730 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"]] 1713 1744 1714 1745 #Can add tide so plots are all on the same tide height … … 1718 1749 max_ele=-100000 1719 1750 min_ele=100000 1720 elevation = [float(x) for x in attribute_dic[" Elevation"]]1751 elevation = [float(x) for x in attribute_dic["elevation"]] 1721 1752 1722 1753 min_ele, max_ele = get_min_max_values(elevation, … … 1743 1774 #print k,l, legend_list_dic[j] 1744 1775 1745 if " Stage" in quantities:1776 if "stage" in quantities: 1746 1777 pylab.figure(100+j) 1747 1778 pylab.plot(t, stage, c = cstr[i], linewidth=1) … … 1750 1781 pylab.axis(stage_axis) 1751 1782 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, 1753 1785 base_name+number, 1754 1786 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 1755 1792 pylab.savefig(figname) 1756 if "Speed" in quantities: 1793 1794 if "speed" in quantities: 1757 1795 pylab.figure(200+j) 1758 1796 pylab.plot(t, speed, c = cstr[i], linewidth=1) … … 1761 1799 pylab.axis(speed_axis) 1762 1800 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, 1764 1803 base_name+number, 1765 1804 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 1766 1810 pylab.savefig(figname) 1767 if " Momentum" in quantities:1811 if "xmomentum" in quantities: 1768 1812 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) 1769 1849 pylab.plot(t, momentum, c = cstr[i], linewidth=1) 1770 1850 pylab.xlabel(time_label) … … 1772 1852 pylab.axis(momentum_axis) 1773 1853 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, 1775 1856 base_name+number, 1776 1857 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 1777 1863 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 1778 1901 if verbose: print 'Closing all plots' 1779 1902 pylab.close('all') 1780 1903 del pylab 1781 1904 if verbose: print 'Finished closing plots' 1905 1906 def 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 1782 2193 1783 2194 def get_min_max_values(list=None,min1=100,max1=-100): … … 1848 2259 file.close() 1849 2260 1850 2261 def 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 1851 2281 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.