Changeset 7866


Ignore:
Timestamp:
Jun 22, 2010, 12:04:32 PM (10 years ago)
Author:
hudson
Message:

More swb tests passing. Cleaned up some pylint errors.

Location:
trunk/anuga_core/source/anuga
Files:
1 deleted
18 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/__init__.py

    r7858 r7866  
    9393# Forcing
    9494#-----------------------------
    95 from anuga.shallow_water.forcing import Inflow
     95from anuga.shallow_water.forcing import Inflow, Rainfall
    9696
    9797#-----------------------------
  • trunk/anuga_core/source/anuga/file/test_csv.py

    r7772 r7866  
    220220       
    221221
     222   
     223    def test_csv2polygons(self):
     224        """test loading of a csv polygon file.
     225        """
     226       
     227        path = get_pathname_from_package('anuga.shallow_water')               
     228        testfile = os.path.join(path, 'polygon_values_example.csv')               
     229
     230        polygons, values = load_csv_as_polygons(testfile,
     231                                        value_name='floors')
     232
     233        assert len(polygons) == 7, 'Must have seven polygons'
     234        assert len(values) == 7, 'Must have seven values'
     235           
     236        # Known floor values
     237        floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
     238       
     239        # Known polygon values
     240        known_polys = {'1': [[422681.61,871117.55],
     241                             [422691.02,871117.60],
     242                             [422690.87,871084.23],
     243                             [422649.36,871081.85],
     244                             [422649.36,871080.39],
     245                             [422631.86,871079.50],
     246                             [422631.72,871086.75],
     247                             [422636.75,871087.20],
     248                             [422636.75,871091.50],
     249                             [422649.66,871092.09],
     250                             [422649.83,871084.91],
     251                             [422652.94,871084.90],
     252                             [422652.84,871092.39],
     253                             [422681.83,871093.73],
     254                             [422681.61,871117.55]],
     255                       '2': [[422664.22,870785.46],
     256                             [422672.48,870780.14],
     257                             [422668.17,870772.62],
     258                             [422660.35,870777.17],
     259                             [422664.22,870785.46]],
     260                       '3': [[422661.30,871215.06],
     261                             [422667.50,871215.70],
     262                             [422668.30,871204.86],
     263                             [422662.21,871204.33],
     264                             [422661.30,871215.06]],
     265                       '4': [[422473.44,871191.22],
     266                             [422478.33,871192.26],
     267                             [422479.52,871186.03],
     268                             [422474.78,871185.14],
     269                             [422473.44,871191.22]],
     270                       '5': [[422369.69,871049.29],
     271                             [422378.63,871053.58],
     272                             [422383.91,871044.51],
     273                             [422374.97,871040.32],
     274                             [422369.69,871049.29]],
     275                       '8': [[422730.56,871203.13],
     276                             [422734.10,871204.90],
     277                             [422735.26,871202.18],
     278                             [422731.87,871200.58],
     279                             [422730.56,871203.13]],
     280                       '9': [[422659.85,871213.80],
     281                             [422660.91,871210.97],
     282                             [422655.42,871208.85],
     283                             [422654.36,871211.68],
     284                             [422659.85,871213.80]]
     285                       }       
     286       
     287
     288       
     289               
     290        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
     291            assert id in polygons.keys()
     292            assert id in values.keys()           
     293
     294            assert int(values[id]) == int(floors[id])
     295            assert len(polygons[id]) == len(known_polys[id])
     296            assert num.allclose(polygons[id], known_polys[id])
     297
     298
     299    def test_csv2polygons_with_clipping(self):
     300        """test_csv2polygons with optional clipping
     301        """
     302        #FIXME(Ole): Not Done!!
     303       
     304        path = get_pathname_from_package('anuga.shallow_water')               
     305        testfile = os.path.join(path, 'polygon_values_example.csv')               
     306
     307        polygons, values = load_csv_as_polygons(testfile,
     308                                        value_name='floors',
     309                                        clipping_polygons=None)
     310
     311        assert len(polygons) == 7, 'Must have seven polygons'
     312        assert len(values) == 7, 'Must have seven values'
     313           
     314        # Known floor values
     315        floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
     316       
     317        # Known polygon values
     318        known_polys = {'1': [[422681.61,871117.55],
     319                             [422691.02,871117.60],
     320                             [422690.87,871084.23],
     321                             [422649.36,871081.85],
     322                             [422649.36,871080.39],
     323                             [422631.86,871079.50],
     324                             [422631.72,871086.75],
     325                             [422636.75,871087.20],
     326                             [422636.75,871091.50],
     327                             [422649.66,871092.09],
     328                             [422649.83,871084.91],
     329                             [422652.94,871084.90],
     330                             [422652.84,871092.39],
     331                             [422681.83,871093.73],
     332                             [422681.61,871117.55]],
     333                       '2': [[422664.22,870785.46],
     334                             [422672.48,870780.14],
     335                             [422668.17,870772.62],
     336                             [422660.35,870777.17],
     337                             [422664.22,870785.46]],
     338                       '3': [[422661.30,871215.06],
     339                             [422667.50,871215.70],
     340                             [422668.30,871204.86],
     341                             [422662.21,871204.33],
     342                             [422661.30,871215.06]],
     343                       '4': [[422473.44,871191.22],
     344                             [422478.33,871192.26],
     345                             [422479.52,871186.03],
     346                             [422474.78,871185.14],
     347                             [422473.44,871191.22]],
     348                       '5': [[422369.69,871049.29],
     349                             [422378.63,871053.58],
     350                             [422383.91,871044.51],
     351                             [422374.97,871040.32],
     352                             [422369.69,871049.29]],
     353                       '8': [[422730.56,871203.13],
     354                             [422734.10,871204.90],
     355                             [422735.26,871202.18],
     356                             [422731.87,871200.58],
     357                             [422730.56,871203.13]],
     358                       '9': [[422659.85,871213.80],
     359                             [422660.91,871210.97],
     360                             [422655.42,871208.85],
     361                             [422654.36,871211.68],
     362                             [422659.85,871213.80]]
     363                       }       
     364       
     365
     366       
     367               
     368        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
     369            assert id in polygons.keys()
     370            assert id in values.keys()           
     371
     372            assert int(values[id]) == int(floors[id])
     373            assert len(polygons[id]) == len(known_polys[id])
     374            assert num.allclose(polygons[id], known_polys[id])
     375
     376
     377
     378
     379   
     380    def test_csv2building_polygons(self):
     381        """test_csv2building_polygons
     382        """
     383       
     384        path = get_pathname_from_package('anuga.shallow_water')               
     385        testfile = os.path.join(path, 'polygon_values_example.csv')               
     386
     387        polygons, values = load_csv_as_building_polygons(testfile,
     388                                                 floor_height=3)
     389
     390        assert len(polygons) == 7, 'Must have seven polygons'
     391        assert len(values) == 7, 'Must have seven values'
     392           
     393        # Known floor values
     394        floors = {'1': 6, '2': 0, '3': 3, '4': 6, '5': 0, '8': 3, '9': 3}
     395       
     396               
     397        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
     398            assert id in polygons.keys()
     399            assert id in values.keys()           
     400
     401            assert float(values[id]) == float(floors[id])
     402
    222403
    223404
  • trunk/anuga_core/source/anuga/file/test_sww.py

    r7770 r7866  
    251251                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
    252252
     253
     254    def test_triangulation(self):
     255        #
     256        # 
     257       
     258        filename = tempfile.mktemp("_data_manager.sww")
     259        outfile = NetCDFFile(filename, netcdf_mode_w)
     260        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
     261        volumes = (0,1,2)
     262        elevation = [0,1,2]
     263        new_origin = None
     264        new_origin = Geo_reference(56, 0, 0)
     265        times = [0, 10]
     266        number_of_volumes = len(volumes)
     267        number_of_points = len(points_utm)
     268        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
     269        sww.store_header(outfile, times, number_of_volumes,
     270                         number_of_points, description='fully sick testing',
     271                         verbose=self.verbose,sww_precision=netcdf_float)
     272        sww.store_triangulation(outfile, points_utm, volumes,
     273                                elevation,  new_origin=new_origin,
     274                                verbose=self.verbose)       
     275        outfile.close()
     276        fid = NetCDFFile(filename)
     277
     278        x = fid.variables['x'][:]
     279        y = fid.variables['y'][:]
     280        fid.close()
     281
     282        assert num.allclose(num.array(map(None, x,y)), points_utm)
     283        os.remove(filename)
     284
     285       
     286    def test_triangulationII(self):
     287        #
     288        # 
     289       
     290        filename = tempfile.mktemp("_data_manager.sww")
     291        outfile = NetCDFFile(filename, netcdf_mode_w)
     292        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
     293        volumes = (0,1,2)
     294        elevation = [0,1,2]
     295        new_origin = None
     296        #new_origin = Geo_reference(56, 0, 0)
     297        times = [0, 10]
     298        number_of_volumes = len(volumes)
     299        number_of_points = len(points_utm)
     300        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
     301        sww.store_header(outfile, times, number_of_volumes,
     302                         number_of_points, description='fully sick testing',
     303                         verbose=self.verbose,sww_precision=netcdf_float)
     304        sww.store_triangulation(outfile, points_utm, volumes,
     305                                new_origin=new_origin,
     306                                verbose=self.verbose)
     307        sww.store_static_quantities(outfile, elevation=elevation)                               
     308                               
     309        outfile.close()
     310        fid = NetCDFFile(filename)
     311
     312        x = fid.variables['x'][:]
     313        y = fid.variables['y'][:]
     314        results_georef = Geo_reference()
     315        results_georef.read_NetCDF(fid)
     316        assert results_georef == Geo_reference(DEFAULT_ZONE, 0, 0)
     317        fid.close()
     318
     319        assert num.allclose(num.array(map(None, x,y)), points_utm)
     320        os.remove(filename)
     321
     322       
     323    def test_triangulation_new_origin(self):
     324        #
     325        # 
     326       
     327        filename = tempfile.mktemp("_data_manager.sww")
     328        outfile = NetCDFFile(filename, netcdf_mode_w)
     329        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
     330        volumes = (0,1,2)
     331        elevation = [0,1,2]
     332        new_origin = None
     333        new_origin = Geo_reference(56, 1, 554354)
     334        points_utm = new_origin.change_points_geo_ref(points_utm)
     335        times = [0, 10]
     336        number_of_volumes = len(volumes)
     337        number_of_points = len(points_utm)
     338        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
     339        sww.store_header(outfile, times, number_of_volumes,
     340                         number_of_points, description='fully sick testing',
     341                         verbose=self.verbose,sww_precision=netcdf_float)
     342        sww.store_triangulation(outfile, points_utm, volumes,
     343                                elevation,  new_origin=new_origin,
     344                                verbose=self.verbose)
     345        outfile.close()
     346        fid = NetCDFFile(filename)
     347
     348        x = fid.variables['x'][:]
     349        y = fid.variables['y'][:]
     350        results_georef = Geo_reference()
     351        results_georef.read_NetCDF(fid)
     352        assert results_georef == new_origin
     353        fid.close()
     354
     355        absolute = Geo_reference(56, 0,0)
     356        assert num.allclose(num.array(
     357            absolute.change_points_geo_ref(map(None, x,y),
     358                                           new_origin)),points_utm)
     359       
     360        os.remove(filename)
     361       
     362    def test_triangulation_points_georeference(self):
     363        #
     364        # 
     365       
     366        filename = tempfile.mktemp("_data_manager.sww")
     367        outfile = NetCDFFile(filename, netcdf_mode_w)
     368        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
     369        volumes = (0,1,2)
     370        elevation = [0,1,2]
     371        new_origin = None
     372        points_georeference = Geo_reference(56, 1, 554354)
     373        points_utm = points_georeference.change_points_geo_ref(points_utm)
     374        times = [0, 10]
     375        number_of_volumes = len(volumes)
     376        number_of_points = len(points_utm)
     377        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
     378        sww.store_header(outfile, times, number_of_volumes,
     379                         number_of_points, description='fully sick testing',
     380                         verbose=self.verbose,sww_precision=netcdf_float)
     381        sww.store_triangulation(outfile, points_utm, volumes,
     382                                elevation,  new_origin=new_origin,
     383                                points_georeference=points_georeference,
     384                                verbose=self.verbose)       
     385        outfile.close()
     386        fid = NetCDFFile(filename)
     387
     388        x = fid.variables['x'][:]
     389        y = fid.variables['y'][:]
     390        results_georef = Geo_reference()
     391        results_georef.read_NetCDF(fid)
     392        assert results_georef == points_georeference
     393        fid.close()
     394
     395        assert num.allclose(num.array(map(None, x,y)), points_utm)
     396        os.remove(filename)
     397       
     398    def test_triangulation_2_geo_refs(self):
     399        #
     400        # 
     401       
     402        filename = tempfile.mktemp("_data_manager.sww")
     403        outfile = NetCDFFile(filename, netcdf_mode_w)
     404        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
     405        volumes = (0,1,2)
     406        elevation = [0,1,2]
     407        new_origin = Geo_reference(56, 1, 1)
     408        points_georeference = Geo_reference(56, 0, 0)
     409        points_utm = points_georeference.change_points_geo_ref(points_utm)
     410        times = [0, 10]
     411        number_of_volumes = len(volumes)
     412        number_of_points = len(points_utm)
     413        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
     414        sww.store_header(outfile, times, number_of_volumes,
     415                         number_of_points, description='fully sick testing',
     416                         verbose=self.verbose,sww_precision=netcdf_float)
     417        sww.store_triangulation(outfile, points_utm, volumes,
     418                                elevation,  new_origin=new_origin,
     419                                points_georeference=points_georeference,
     420                                verbose=self.verbose)       
     421        outfile.close()
     422        fid = NetCDFFile(filename)
     423
     424        x = fid.variables['x'][:]
     425        y = fid.variables['y'][:]
     426        results_georef = Geo_reference()
     427        results_georef.read_NetCDF(fid)
     428        assert results_georef == new_origin
     429        fid.close()
     430
     431
     432        absolute = Geo_reference(56, 0,0)
     433        assert num.allclose(num.array(
     434            absolute.change_points_geo_ref(map(None, x,y),
     435                                           new_origin)),points_utm)
     436        os.remove(filename)
     437
    253438#################################################################################
    254439
  • trunk/anuga_core/source/anuga/file_conversion/file_conversion.py

    r7814 r7866  
    9191# @param quantity_names
    9292# @param time_as_seconds
    93 def timefile2netcdf(file_text, quantity_names=None, time_as_seconds=False):
     93def timefile2netcdf(file_text, file_out = None, quantity_names=None, \
     94                                time_as_seconds=False):
    9495    """Template for converting typical text files with time series to
    9596    NetCDF tms file.
     
    111112    will provide a time dependent function f(t) with three attributes
    112113
    113     filename is assumed to be the rootname with extenisons .txt and .sww
     114    filename is assumed to be the rootname with extensions .txt/.tms and .sww
    114115    """
    115116
     
    121122        raise IOError('Input file %s should be of type .txt.' % file_text)
    122123
    123     filename = file_text[:-4]
     124    if file_out == None:
     125        file_out = file_text[:-4] + '.tms'
     126
    124127    fid = open(file_text)
    125128    line = fid.readline()
     
    181184            Q[i, j] = float(value)
    182185
    183     msg = 'File %s must list time as a monotonuosly ' % filename
     186    msg = 'File %s must list time as a monotonuosly ' % file_text
    184187    msg += 'increasing sequence'
    185188    assert num.alltrue(T[1:] - T[:-1] > 0), msg
     
    188191    from Scientific.IO.NetCDF import NetCDFFile
    189192
    190     fid = NetCDFFile(filename + '.tms', netcdf_mode_w)
     193    fid = NetCDFFile(file_out, netcdf_mode_w)
    191194
    192195    fid.institution = 'Geoscience Australia'
  • trunk/anuga_core/source/anuga/file_conversion/test_sww2dem.py

    r7841 r7866  
    15201520        os.remove(ascfile)
    15211521        os.remove(swwfile2)
     1522
     1523
     1524    def test_export_grid(self):
     1525        """
     1526        test_export_grid(self):
     1527        Test that sww information can be converted correctly to asc/prj
     1528        format readable by e.g. ArcView
     1529        """
     1530
     1531        import time, os
     1532        from Scientific.IO.NetCDF import NetCDFFile
     1533
     1534        try:
     1535            os.remove('teg*.sww')
     1536        except:
     1537            pass
     1538
     1539        #Setup
     1540        self.domain.set_name('teg')
     1541
     1542        prjfile = self.domain.get_name() + '_elevation.prj'
     1543        ascfile = self.domain.get_name() + '_elevation.asc'
     1544        swwfile = self.domain.get_name() + '.sww'
     1545
     1546        self.domain.set_datadir('.')
     1547        self.domain.smooth = True
     1548        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     1549        self.domain.set_quantity('stage', 1.0)
     1550
     1551        self.domain.geo_reference = Geo_reference(56,308500,6189000)
     1552
     1553        sww = SWW_file(self.domain)
     1554        sww.store_connectivity()
     1555        sww.store_timestep()
     1556        self.domain.evolve_to_end(finaltime = 0.01)
     1557        sww.store_timestep()
     1558
     1559        cellsize = 0.25
     1560        #Check contents
     1561        #Get NetCDF
     1562
     1563        fid = NetCDFFile(sww.filename, netcdf_mode_r)
     1564
     1565        # Get the variables
     1566        x = fid.variables['x'][:]
     1567        y = fid.variables['y'][:]
     1568        z = fid.variables['elevation'][:]
     1569        time = fid.variables['time'][:]
     1570        stage = fid.variables['stage'][:]
     1571
     1572        fid.close()
     1573
     1574        #Export to ascii/prj files
     1575        sww2dem_batch(self.domain.get_name(),
     1576                quantities = 'elevation',
     1577                cellsize = cellsize,
     1578                verbose = self.verbose,
     1579                format = 'asc')
     1580
     1581        #Check asc file
     1582        ascid = open(ascfile)
     1583        lines = ascid.readlines()
     1584        ascid.close()
     1585
     1586        L = lines[2].strip().split()
     1587        assert L[0].strip().lower() == 'xllcorner'
     1588        assert num.allclose(float(L[1].strip().lower()), 308500)
     1589
     1590        L = lines[3].strip().split()
     1591        assert L[0].strip().lower() == 'yllcorner'
     1592        assert num.allclose(float(L[1].strip().lower()), 6189000)
     1593
     1594        #Check grid values
     1595        for j in range(5):
     1596            L = lines[6+j].strip().split()
     1597            y = (4-j) * cellsize
     1598            for i in range(5):
     1599                assert num.allclose(float(L[i]), -i*cellsize - y)
     1600               
     1601        #Cleanup
     1602        os.remove(prjfile)
     1603        os.remove(ascfile)
     1604        os.remove(swwfile)
     1605
     1606    def test_export_gridII(self):
     1607        """
     1608        test_export_gridII(self):
     1609        Test that sww information can be converted correctly to asc/prj
     1610        format readable by e.g. ArcView
     1611        """
     1612
     1613        import time, os
     1614        from Scientific.IO.NetCDF import NetCDFFile
     1615
     1616        try:
     1617            os.remove('teg*.sww')
     1618        except:
     1619            pass
     1620
     1621        #Setup
     1622        self.domain.set_name('tegII')
     1623
     1624        swwfile = self.domain.get_name() + '.sww'
     1625
     1626        self.domain.set_datadir('.')
     1627        self.domain.smooth = True
     1628        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     1629        self.domain.set_quantity('stage', 1.0)
     1630
     1631        self.domain.geo_reference = Geo_reference(56,308500,6189000)
     1632
     1633        sww = SWW_file(self.domain)
     1634        sww.store_connectivity()
     1635        sww.store_timestep()
     1636        self.domain.evolve_to_end(finaltime = 0.01)
     1637        sww.store_timestep()
     1638
     1639        cellsize = 0.25
     1640        #Check contents
     1641        #Get NetCDF
     1642
     1643        fid = NetCDFFile(sww.filename, netcdf_mode_r)
     1644
     1645        # Get the variables
     1646        x = fid.variables['x'][:]
     1647        y = fid.variables['y'][:]
     1648        z = fid.variables['elevation'][:]
     1649        time = fid.variables['time'][:]
     1650        stage = fid.variables['stage'][:]
     1651        xmomentum = fid.variables['xmomentum'][:]
     1652        ymomentum = fid.variables['ymomentum'][:]       
     1653
     1654        #print 'stage', stage
     1655        #print 'xmom', xmomentum
     1656        #print 'ymom', ymomentum       
     1657
     1658        fid.close()
     1659
     1660        #Export to ascii/prj files
     1661        if True:
     1662            sww2dem_batch(self.domain.get_name(),
     1663                        quantities = ['elevation', 'depth'],
     1664                        cellsize = cellsize,
     1665                        verbose = self.verbose,
     1666                        format = 'asc')
     1667
     1668        else:
     1669            sww2dem_batch(self.domain.get_name(),
     1670                quantities = ['depth'],
     1671                cellsize = cellsize,
     1672                verbose = self.verbose,
     1673                format = 'asc')
     1674
     1675
     1676            export_grid(self.domain.get_name(),
     1677                quantities = ['elevation'],
     1678                cellsize = cellsize,
     1679                verbose = self.verbose,
     1680                format = 'asc')
     1681
     1682        prjfile = self.domain.get_name() + '_elevation.prj'
     1683        ascfile = self.domain.get_name() + '_elevation.asc'
     1684       
     1685        #Check asc file
     1686        ascid = open(ascfile)
     1687        lines = ascid.readlines()
     1688        ascid.close()
     1689
     1690        L = lines[2].strip().split()
     1691        assert L[0].strip().lower() == 'xllcorner'
     1692        assert num.allclose(float(L[1].strip().lower()), 308500)
     1693
     1694        L = lines[3].strip().split()
     1695        assert L[0].strip().lower() == 'yllcorner'
     1696        assert num.allclose(float(L[1].strip().lower()), 6189000)
     1697
     1698        #print "ascfile", ascfile
     1699        #Check grid values
     1700        for j in range(5):
     1701            L = lines[6+j].strip().split()
     1702            y = (4-j) * cellsize
     1703            for i in range(5):
     1704                #print " -i*cellsize - y",  -i*cellsize - y
     1705                #print "float(L[i])", float(L[i])
     1706                assert num.allclose(float(L[i]), -i*cellsize - y)
     1707
     1708        #Cleanup
     1709        os.remove(prjfile)
     1710        os.remove(ascfile)
     1711       
     1712        #Check asc file
     1713        ascfile = self.domain.get_name() + '_depth.asc'
     1714        prjfile = self.domain.get_name() + '_depth.prj'
     1715        ascid = open(ascfile)
     1716        lines = ascid.readlines()
     1717        ascid.close()
     1718
     1719        L = lines[2].strip().split()
     1720        assert L[0].strip().lower() == 'xllcorner'
     1721        assert num.allclose(float(L[1].strip().lower()), 308500)
     1722
     1723        L = lines[3].strip().split()
     1724        assert L[0].strip().lower() == 'yllcorner'
     1725        assert num.allclose(float(L[1].strip().lower()), 6189000)
     1726
     1727        #Check grid values
     1728        for j in range(5):
     1729            L = lines[6+j].strip().split()
     1730            y = (4-j) * cellsize
     1731            for i in range(5):
     1732                #print " -i*cellsize - y",  -i*cellsize - y
     1733                #print "float(L[i])", float(L[i])               
     1734                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
     1735
     1736        #Cleanup
     1737        os.remove(prjfile)
     1738        os.remove(ascfile)
     1739        os.remove(swwfile)
     1740
     1741
     1742    def test_export_gridIII(self):
     1743        """
     1744        test_export_gridIII
     1745        Test that sww information can be converted correctly to asc/prj
     1746        format readable by e.g. ArcView
     1747        """
     1748
     1749        import time, os
     1750        from Scientific.IO.NetCDF import NetCDFFile
     1751
     1752        try:
     1753            os.remove('teg*.sww')
     1754        except:
     1755            pass
     1756
     1757        #Setup
     1758       
     1759        self.domain.set_name('tegIII')
     1760
     1761        swwfile = self.domain.get_name() + '.sww'
     1762
     1763        self.domain.set_datadir('.')
     1764        self.domain.format = 'sww'
     1765        self.domain.smooth = True
     1766        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     1767        self.domain.set_quantity('stage', 1.0)
     1768
     1769        self.domain.geo_reference = Geo_reference(56,308500,6189000)
     1770       
     1771        sww = SWW_file(self.domain)
     1772        sww.store_connectivity()
     1773        sww.store_timestep() #'stage')
     1774        self.domain.evolve_to_end(finaltime = 0.01)
     1775        sww.store_timestep() #'stage')
     1776
     1777        cellsize = 0.25
     1778        #Check contents
     1779        #Get NetCDF
     1780
     1781        fid = NetCDFFile(sww.filename, netcdf_mode_r)
     1782
     1783        # Get the variables
     1784        x = fid.variables['x'][:]
     1785        y = fid.variables['y'][:]
     1786        z = fid.variables['elevation'][:]
     1787        time = fid.variables['time'][:]
     1788        stage = fid.variables['stage'][:]
     1789
     1790        fid.close()
     1791
     1792        #Export to ascii/prj files
     1793        extra_name_out = 'yeah'
     1794        if True:
     1795            sww2dem_batch(self.domain.get_name(),
     1796                        quantities = ['elevation', 'depth'],
     1797                        extra_name_out = extra_name_out,
     1798                        cellsize = cellsize,
     1799                        verbose = self.verbose,
     1800                        format = 'asc')
     1801
     1802        else:
     1803            sww2dem_batch(self.domain.get_name(),
     1804                quantities = ['depth'],
     1805                cellsize = cellsize,
     1806                verbose = self.verbose,
     1807                format = 'asc')
     1808
     1809
     1810            sww2dem_batch(self.domain.get_name(),
     1811                quantities = ['elevation'],
     1812                cellsize = cellsize,
     1813                verbose = self.verbose,
     1814                format = 'asc')
     1815
     1816        prjfile = self.domain.get_name() + '_elevation_yeah.prj'
     1817        ascfile = self.domain.get_name() + '_elevation_yeah.asc'
     1818       
     1819        #Check asc file
     1820        ascid = open(ascfile)
     1821        lines = ascid.readlines()
     1822        ascid.close()
     1823
     1824        L = lines[2].strip().split()
     1825        assert L[0].strip().lower() == 'xllcorner'
     1826        assert num.allclose(float(L[1].strip().lower()), 308500)
     1827
     1828        L = lines[3].strip().split()
     1829        assert L[0].strip().lower() == 'yllcorner'
     1830        assert num.allclose(float(L[1].strip().lower()), 6189000)
     1831
     1832        #print "ascfile", ascfile
     1833        #Check grid values
     1834        for j in range(5):
     1835            L = lines[6+j].strip().split()
     1836            y = (4-j) * cellsize
     1837            for i in range(5):
     1838                #print " -i*cellsize - y",  -i*cellsize - y
     1839                #print "float(L[i])", float(L[i])
     1840                assert num.allclose(float(L[i]), -i*cellsize - y)
     1841               
     1842        #Cleanup
     1843        os.remove(prjfile)
     1844        os.remove(ascfile)
     1845       
     1846        #Check asc file
     1847        ascfile = self.domain.get_name() + '_depth_yeah.asc'
     1848        prjfile = self.domain.get_name() + '_depth_yeah.prj'
     1849        ascid = open(ascfile)
     1850        lines = ascid.readlines()
     1851        ascid.close()
     1852
     1853        L = lines[2].strip().split()
     1854        assert L[0].strip().lower() == 'xllcorner'
     1855        assert num.allclose(float(L[1].strip().lower()), 308500)
     1856
     1857        L = lines[3].strip().split()
     1858        assert L[0].strip().lower() == 'yllcorner'
     1859        assert num.allclose(float(L[1].strip().lower()), 6189000)
     1860
     1861        #Check grid values
     1862        for j in range(5):
     1863            L = lines[6+j].strip().split()
     1864            y = (4-j) * cellsize
     1865            for i in range(5):
     1866                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
     1867
     1868        #Cleanup
     1869        os.remove(prjfile)
     1870        os.remove(ascfile)
     1871        os.remove(swwfile)
     1872
     1873    def test_export_grid_bad(self):
     1874        """Test that sww information can be converted correctly to asc/prj
     1875        format readable by e.g. ArcView
     1876        """
     1877
     1878        try:
     1879            sww2dem_batch('a_small_round-egg',
     1880                        quantities = ['elevation', 'depth'],
     1881                        cellsize = 99,
     1882                        verbose = self.verbose,
     1883                        format = 'asc')
     1884        except IOError:
     1885            pass
     1886        else:
     1887            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
    15221888       
    15231889
  • trunk/anuga_core/source/anuga/file_conversion/test_urs2sts.py

    r7847 r7866  
    19581958        os.remove(meshname)
    19591959       
    1960 
     1960       
     1961       
     1962    def test_file_boundary_sts_time_limit(self):
     1963        """test_file_boundary_stsIV_sinewave_ordering(self):
     1964        Read correct points from ordering file and apply sts to boundary
     1965        This one uses a sine wave and compares to time boundary
     1966       
     1967        This one tests that times used can be limited by upper limit
     1968        """
     1969       
     1970        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
     1971        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
     1972        tide = 0.35
     1973        time_step_count = 50
     1974        time_step = 0.1
     1975        times_ref = num.arange(0, time_step_count*time_step, time_step)
     1976       
     1977        n=len(lat_long_points)
     1978        first_tstep=num.ones(n,num.int)
     1979        last_tstep=(time_step_count)*num.ones(n,num.int)
     1980       
     1981        gauge_depth=20*num.ones(n,num.float)
     1982       
     1983        ha1=num.ones((n,time_step_count),num.float)
     1984        ua1=3.*num.ones((n,time_step_count),num.float)
     1985        va1=2.*num.ones((n,time_step_count),num.float)
     1986        for i in range(n):
     1987            ha1[i]=num.sin(times_ref)
     1988       
     1989       
     1990        base_name, files = self.write_mux2(lat_long_points,
     1991                                           time_step_count, time_step,
     1992                                           first_tstep, last_tstep,
     1993                                           depth=gauge_depth,
     1994                                           ha=ha1,
     1995                                           ua=ua1,
     1996                                           va=va1)
     1997
     1998        # Write order file
     1999        file_handle, order_base_name = tempfile.mkstemp("")
     2000        os.close(file_handle)
     2001        os.remove(order_base_name)
     2002        d=","
     2003        order_file=order_base_name+'order.txt'
     2004        fid=open(order_file,'w')
     2005       
     2006        # Write Header
     2007        header='index, longitude, latitude\n'
     2008        fid.write(header)
     2009        indices=[3,0,1]
     2010        for i in indices:
     2011            line=str(i)+d+str(lat_long_points[i][1])+d+\
     2012                str(lat_long_points[i][0])+"\n"
     2013            fid.write(line)
     2014        fid.close()
     2015
     2016        sts_file=base_name
     2017        urs2sts(base_name, basename_out=sts_file,
     2018                ordering_filename=order_file,
     2019                mean_stage=tide,
     2020                verbose=False)
     2021        self.delete_mux(files)
     2022       
     2023       
     2024       
     2025        # Now read the sts file and check that values have been stored correctly.
     2026        fid = NetCDFFile(sts_file + '.sts')
     2027
     2028        # Check the time vector
     2029        times = fid.variables['time'][:]
     2030        starttime = fid.starttime[0]
     2031        #print times
     2032        #print starttime
     2033
     2034        # Check sts quantities
     2035        stage = fid.variables['stage'][:]
     2036        xmomentum = fid.variables['xmomentum'][:]
     2037        ymomentum = fid.variables['ymomentum'][:]
     2038        elevation = fid.variables['elevation'][:]
     2039
     2040       
     2041
     2042        # Create beginnings of boundary polygon based on sts_boundary
     2043        boundary_polygon = create_sts_boundary(base_name)
     2044       
     2045        os.remove(order_file)
     2046
     2047        # Append the remaining part of the boundary polygon to be defined by
     2048        # the user
     2049        bounding_polygon_utm=[]
     2050        for point in bounding_polygon:
     2051            zone,easting,northing=redfearn(point[0],point[1])
     2052            bounding_polygon_utm.append([easting,northing])
     2053
     2054        boundary_polygon.append(bounding_polygon_utm[3])
     2055        boundary_polygon.append(bounding_polygon_utm[4])
     2056
     2057        #print 'boundary_polygon', boundary_polygon
     2058       
     2059
     2060        assert num.allclose(bounding_polygon_utm,boundary_polygon)
     2061
     2062
     2063        extent_res=1000000
     2064        meshname = 'urs_test_mesh' + '.tsh'
     2065        interior_regions=None
     2066        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
     2067       
     2068        # have to change boundary tags from last example because now bounding
     2069        # polygon starts in different place.
     2070        create_mesh_from_regions(boundary_polygon,
     2071                                 boundary_tags=boundary_tags,
     2072                                 maximum_triangle_area=extent_res,
     2073                                 filename=meshname,
     2074                                 interior_regions=interior_regions,
     2075                                 verbose=False)
     2076       
     2077        domain_fbound = Domain(meshname)
     2078        domain_fbound.set_quantity('stage', tide)
     2079       
     2080       
     2081        Bf = File_boundary(sts_file+'.sts',
     2082                           domain_fbound,
     2083                           boundary_polygon=boundary_polygon)
     2084        time_vec = Bf.F.get_time()
     2085        assert num.allclose(Bf.F.starttime, starttime)
     2086        assert num.allclose(time_vec, times_ref)                                   
     2087       
     2088        for time_limit in [0.1, 0.2, 0.5, 1.0, 2.2, 3.0, 4.3, 6.0, 10.0]:
     2089            Bf = File_boundary(sts_file+'.sts',
     2090                               domain_fbound,
     2091                               time_limit=time_limit+starttime,
     2092                               boundary_polygon=boundary_polygon)
     2093       
     2094            time_vec = Bf.F.get_time()
     2095            assert num.allclose(Bf.F.starttime, starttime)           
     2096            assert num.alltrue(time_vec < time_limit)
     2097           
     2098           
     2099        try:   
     2100            Bf = File_boundary(sts_file+'.sts',
     2101                               domain_fbound,
     2102                               time_limit=-1+starttime,
     2103                               boundary_polygon=boundary_polygon)           
     2104            time_vec = Bf.F.get_time()   
     2105            print time_vec   
     2106        except AssertionError:
     2107            pass
     2108        else:
     2109            raise Exception, 'Should have raised Exception here'
    19612110
    19622111#-------------------------------------------------------------
  • trunk/anuga_core/source/anuga/geometry/quad.py

    r7858 r7866  
    1313"""
    1414
    15 from anuga.utilities.treenode import TreeNode
    1615import anuga.utilities.log as log
    1716
    1817           
    19 class Cell(TreeNode):
     18class Cell():
    2019    """class Cell
    2120
     
    2928            parent is the node above this one, or None if it is root.
    3029        """
    31  
    32         # Initialise base classes
    33         TreeNode.__init__(self, name)
    3430   
    3531        self.extents = extents
     
    4743            ret_str += ', children: %d' % (len(self.children))
    4844        return ret_str
    49 
    50    
    51 
    52     def clear(self):
    53         """ Remove all leaves from this node.
    54         """
    55         self.Prune()   # TreeNode method
    56 
    57 
    58     def clear_leaf_node(self):
    59         """ Clears storage in leaf node.
    60             Called from Treenode.
    61             Must exist.   
    62         """
    63         self.leaves = []
    64    
    65    
    66     def clear_internal_node(self):
    67         """Called from Treenode.   
    68     Must exist.
    69     """
    70         self.leaves = []
    7145
    7246
  • trunk/anuga_core/source/anuga/geometry/test_geometry.py

    r7720 r7866  
     1""" Test for the geometry classes.
     2
     3    Pylint quality rating as of June 2010: 8.51/10.
     4"""
     5
    16import unittest
    2 import numpy as num
    37
    48from aabb import AABB
    59from quad import Cell
    610
    7 import types, sys
     11import types
    812
    913#-------------------------------------------------------------
    1014
    1115class Test_Geometry(unittest.TestCase):
    12 
     16    """ Test geometry classes. """
    1317    def setUp(self):
     18        """ Generic set up for geometry tests. """
    1419        pass
    1520
    1621    def tearDown(self):
     22        """ Generic shut down for geometry tests. """
    1723        pass
    1824
    19     def test_AABB_contains(self):
     25    def test_aabb_contains(self):
     26        """ Test if point is correctly classified as falling inside or
     27            outside of bounding box. """
    2028        box = AABB(1, 21, 1, 11)
    2129        assert box.contains([10, 5])
     
    2836        assert not box.contains([50, 6])       
    2937       
    30     def test_AABB_split_vert(self):
     38    def test_aabb_split_vert(self):
     39        """ Test that a bounding box can be split correctly along an axis.
     40        """
    3141        parent = AABB(1, 21, 1, 11)
    3242       
     
    4353        self.assertEqual(child2.ymax, 11)   
    4454
    45     def test_AABB_split_horiz(self):
     55    def test_aabb_split_horiz(self):
     56        """ Test that a bounding box will be split along the horizontal axis
     57        correctly. """
    4658        parent = AABB(1, 11, 1, 41)
    4759       
     
    5971       
    6072    def test_add_data(self):
    61         cell = Cell(AABB(0,10, 0,5), None)
    62         cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), 222),  \
     73        """ Test add and retrieve arbitrary data from tree structure. """
     74        cell = Cell(AABB(0, 10, 0, 5), None)
     75        cell.insert([(AABB(1, 3, 1, 3), 111), (AABB(8, 9, 1, 2), 222),  \
    6376                     (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)])
    6477
    6578        result = cell.retrieve()
    66         assert type(result) in [types.ListType,types.TupleType],\
     79        assert type(result) in [types.ListType,types.TupleType], \
    6780                            'should be a list'
    6881
    69         self.assertEqual(len(result),4)
     82        self.assertEqual(len(result), 4)
    7083       
    7184    def test_search(self):
     85        """ Test search tree for an intersection. """
    7286        test_tag = 222
    73         cell = Cell(AABB(0,10, 0,5), None)
    74         cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), test_tag),  \
     87        cell = Cell(AABB(0, 10, 0,5), None)
     88        cell.insert([(AABB(1, 3, 1, 3), 111), (AABB(8, 9, 1, 2), test_tag),  \
    7589                     (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)])
    7690
    7791        result = cell.search([8.5, 1.5])
    78         assert type(result) in [types.ListType,types.TupleType],\
     92        assert type(result) in [types.ListType, types.TupleType], \
    7993                            'should be a list'
    8094        assert(len(result) == 1)
    81         data, node = result[0]
     95        data, _ = result[0]
    8296        self.assertEqual(data, test_tag, 'only 1 point should intersect')
    8397
    8498    def test_get_siblings(self):
    8599        """ Make sure children know their parent. """
    86         cell = Cell(AABB(0,10, 0,5), None)
    87         cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), 222)])
     100        cell = Cell(AABB(0, 10, 0, 5), None)
     101        cell.insert([(AABB(1, 3, 1, 3), 111), (AABB(8, 9, 1, 2), 222)])
    88102
    89103        assert len(cell.children) == 2
     
    92106        assert cell.children[1].parent == cell
    93107
    94     def test_clear_1(self):
    95         cell = Cell(AABB(0,10, 0,5), None)   
    96         cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), 222),  \
    97                      (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)])
    98                      
    99         assert len(cell.retrieve()) == 4
    100         cell.clear()
    101 
    102         assert len(cell.retrieve()) == 0
    103108
    104109################################################################################
  • trunk/anuga_core/source/anuga/geometry/test_polygon.py

    r7858 r7866  
    11#!/usr/bin/env python
     2
     3""" Test suite to test polygon functionality. """
    24
    35import unittest
    46import numpy as num
    5 from math import sqrt, pi
    67from anuga.utilities.numerical_tools import ensure_numeric
    78from anuga.utilities.system_tools import get_pathname_from_package
    89
    9 from polygon import *
    10 from polygon import _poly_xy
     10from polygon import _poly_xy, separate_points_by_polygon, \
     11                    populate_polygon, polygon_area, is_inside_polygon, \
     12                    read_polygon, point_on_line, point_in_polygon, \
     13                    plot_polygons, outside_polygon, is_outside_polygon, \
     14                    intersection, is_complex, is_inside_triangle, \
     15                    interpolate_polyline, inside_polygon, in_and_outside_polygon
     16                   
    1117from polygon_function import Polygon_function
    1218from anuga.coordinate_transforms.geo_reference import Geo_reference
     
    401407
    402408        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    403         points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     409        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], \
     410                        [0.5, -0.5]]
    404411        res, count = separate_points_by_polygon( points, polygon )
    405412
     
    14011408
    14021409    def test_no_intersection(self):
     1410        """ Test 2 non-touching lines don't intersect. """
    14031411        line0 = [[-1,1], [1,1]]
    14041412        line1 = [[0,-1], [0,0]]
     
    14691477
    14701478        # Overlap
    1471         line0 = [[0,5], [7,19]]
    1472         line1 = [[1,7], [10,25]]
     1479        line0 = [[0, 5], [7, 19]]
     1480        line1 = [[1, 7], [10, 25]]
    14731481        status, value = intersection(line0, line1)
    14741482        assert status == 2
     
    15121520
    15131521
    1514     def NOtest_inside_polygon_main(self):
    1515         # FIXME (Ole): Why is this disabled?
    1516         #print "inside",inside
    1517         #print "outside",outside
    1518 
    1519         assert not inside_polygon((0.5, 1.5), polygon)
    1520         assert not inside_polygon((0.5, -0.5), polygon)
    1521         assert not inside_polygon((-0.5, 0.5), polygon)
    1522         assert not inside_polygon((1.5, 0.5), polygon)
    1523 
    1524         # Try point on borders
    1525         assert inside_polygon((1., 0.5), polygon, closed=True)
    1526         assert inside_polygon((0.5, 1), polygon, closed=True)
    1527         assert inside_polygon((0., 0.5), polygon, closed=True)
    1528         assert inside_polygon((0.5, 0.), polygon, closed=True)
    1529 
    1530         assert not inside_polygon((0.5, 1), polygon, closed=False)
    1531         assert not inside_polygon((0., 0.5), polygon, closed=False)
    1532         assert not inside_polygon((0.5, 0.), polygon, closed=False)
    1533         assert not inside_polygon((1., 0.5), polygon, closed=False)
    1534 
     1522    def test_inside_polygon_main(self):
    15351523        # From real example (that failed)
    15361524        polygon = [[20,20], [40,20], [40,40], [20,40]]
     
    15461534
    15471535    def test_polygon_area(self):
     1536        """ Test getting the area of a polygon. """
    15481537        # Simplest case: Polygon is the unit square
    15491538        polygon = [[0,0], [1,0], [1,1], [0,1]]
     
    16111600        # Another case
    16121601        polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
    1613         v = plot_polygons([polygon2,polygon3], figname='test2')
     1602        plot_polygons([polygon2, polygon3], figname='test2')
    16141603
    16151604        for file in ['test1.png', 'test2.png']:
     
    16191608
    16201609    def test_inside_polygon_geospatial(self):
     1610        """ Test geospatial coords inside poly. """
    16211611        #Simplest case: Polygon is the unit square
    1622         polygon_absolute = [[0,0], [1,0], [1,1], [0,1]]
     1612        polygon_absolute = [[0, 0], [1, 0], [1, 1], [0, 1]]
    16231613        poly_geo_ref = Geo_reference(57, 100, 100)
    16241614        polygon = poly_geo_ref.change_points_geo_ref(polygon_absolute)
     
    16391629        assert is_inside_polygon(points_absolute, polygon_absolute)
    16401630
    1641     def NOtest_decimate_polygon(self):
     1631    def test_decimate_polygon(self):
     1632        from polygon import decimate_polygon
    16421633        polygon = [[0,0], [10,10], [15,5], [20, 10],
    16431634                   [25,0], [30,10], [40,-10], [35, -5]]
     
    16451636        dpoly = decimate_polygon(polygon, factor=2)
    16461637
    1647         print dpoly
    1648 
    16491638        assert len(dpoly)*2==len(polygon)
    16501639
    1651         minx = maxx = polygon[0][0]
    1652         miny = maxy = polygon[0][1]
    1653         for point in polygon[1:]:
    1654             x, y = point
    1655 
    1656             if x < minx: minx = x
    1657             if x > maxx: maxx = x
    1658             if y < miny: miny = y
    1659             if y > maxy: maxy = y
    1660 
    1661         assert [minx, miny] in polygon
    1662         print minx, maxy
    1663         assert [minx, maxy] in polygon
    1664         assert [maxx, miny] in polygon
    1665         assert [maxx, maxy] in polygon
    16661640
    16671641    def test_interpolate_polyline(self):
     
    17561730
    17571731    def test_is_inside_triangle_more(self):
    1758 
     1732        """ Test if points inside triangles are detected correctly. """
    17591733        res = is_inside_triangle([0.5, 0.5], [[ 0.5,  0. ],
    17601734                                              [ 0.5,  0.5],
     
    18021776       
    18031777    def test_is_polygon_complex(self):
    1804         """Test a concave and a complex poly with is_complex, to make sure it can detect
    1805            self-intersection.
     1778        """ Test a concave and a complex poly with is_complex, to make
     1779            sure it can detect self-intersection.
    18061780        """
    18071781        concave_poly = [[0, 0], [10, 0], [5, 5], [10, 10], [0, 10]]
    1808         complex_poly = [[0, 0], [10, 0], [5, 5], [4, 15], [5, 7], [10, 10], [0, 10]]
     1782        complex_poly = [[0, 0], [10, 0], [5, 5], [4, 15], [5, 7], [10, 10], \
     1783                            [0, 10]]
    18091784
    18101785        assert not is_complex(concave_poly)
     
    18121787
    18131788    def test_is_polygon_complex2(self):
    1814         """Test a concave and a complex poly with is_complex, to make sure it can detect
    1815            self-intersection. This test uses more complicated polygons.
     1789        """ Test a concave and a complex poly with is_complex, to make sure it
     1790            can detect self-intersection. This test uses more complicated
     1791            polygons.
    18161792        """   
    1817         concave_poly = [[0, 0], [10, 0], [11,0], [5, 5], [7,6], [10, 10], [1,5], [0, 10]]
    1818         complex_poly = [[0, 0], [12,12], [10, 0], [5, 5], [3,18], [4, 15], [5, 7], [10, 10], [0, 10], [16, 12]]       
     1793        concave_poly = [[0, 0], [10, 0], [11,0], [5, 5], [7, 6], [10, 10], \
     1794                        [1, 5], [0, 10]]
     1795        complex_poly = [[0, 0], [12, 12], [10, 0], [5, 5], [3,18], [4, 15], \
     1796                        [5, 7], [10, 10], [0, 10], [16, 12]]       
    18191797
    18201798        assert not is_complex(concave_poly)
  • trunk/anuga_core/source/anuga/pmesh/mesh_quadtree.py

    r7810 r7866  
    1818    from polygon_ext import _is_inside_triangle       
    1919else:
    20     msg = 'C implementations could not be accessed by %s.\n ' % __file__
    21     msg += 'Make sure compile_all.py has been run as described in '
    22     msg += 'the ANUGA installation guide.'
    23     raise Exception, msg
     20    MESSAGE = 'C implementations could not be accessed by %s.\n ' % __file__
     21    MESSAGE += 'Make sure compile_all.py has been run as described in '
     22    MESSAGE += 'the ANUGA installation guide.'
     23    raise Exception(MESSAGE)
    2424
    2525import numpy as num
     
    5050        extents.grow(1.001) # To avoid round off error
    5151        Cell.__init__(self, extents, None)  # root has no parent
    52        
     52
     53        self.last_triangle = None       
    5354        N = len(mesh)
    5455        self.mesh = mesh
  • trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py

    r7841 r7866  
    503503        sww.store_timestep()
    504504
    505 
    506505        # Check contents
    507506        # Get NetCDF
     
    530529        assert num.allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
    531530                                   Q0[5] + Q2[6] + Q1[7])/6)
    532 
    533 
    534         fid.close()
    535 
    536         #Cleanup
    537         os.remove(sww.filename)
    538 
    539     def no_test_sww_variable3(self):
    540         """Test that sww information can be written correctly
    541         multiple timesteps using a different reduction operator (min)
    542         """
    543 
    544         # Different reduction in sww files has been made obsolete.
    545        
    546         import time, os
    547         from Scientific.IO.NetCDF import NetCDFFile
    548 
    549         self.domain.set_name('datatest' + str(id(self)))
    550         self.domain.format = 'sww'
    551         self.domain.smooth = True
    552         self.domain.reduction = min
    553 
    554         sww = SWW_file(self.domain)
    555         sww.store_connectivity()
    556         sww.store_timestep()
    557         #self.domain.tight_slope_limiters = 1
    558         self.domain.evolve_to_end(finaltime = 0.01)
    559         sww.store_timestep()
    560 
    561 
    562         #Check contents
    563         #Get NetCDF
    564         fid = NetCDFFile(sww.filename, netcdf_mode_r)
    565 
    566         # Get the variables
    567         x = fid.variables['x']
    568         y = fid.variables['y']
    569         z = fid.variables['elevation']
    570         time = fid.variables['time']
    571         stage = fid.variables['stage']
    572 
    573         #Check values
    574         Q = self.domain.quantities['stage']
    575         Q0 = Q.vertex_values[:,0]
    576         Q1 = Q.vertex_values[:,1]
    577         Q2 = Q.vertex_values[:,2]
    578 
    579         A = stage[1,:]
    580         assert num.allclose(A[0], min(Q2[0], Q1[1]))
    581         assert num.allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
    582         assert num.allclose(A[2], Q0[3])
    583         assert num.allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
    584 
    585         #Center point
    586         assert num.allclose(A[4], min(Q1[0], Q2[1], Q0[2],
    587                                       Q0[5], Q2[6], Q1[7]))
    588531
    589532
     
    696639
    697640
    698     def Not_a_test_sww_DSG(self):
    699         """Not a test, rather a look at the sww format
    700         """
    701 
    702         import time, os
    703         from Scientific.IO.NetCDF import NetCDFFile
    704 
    705         self.domain.set_name('datatest' + str(id(self)))
    706         self.domain.format = 'sww'
    707         self.domain.smooth = True
    708         self.domain.reduction = mean
    709 
    710         sww = SWW_file(self.domain)
    711         sww.store_connectivity()
    712         sww.store_timestep()
    713 
    714         #Check contents
    715         #Get NetCDF
    716         fid = NetCDFFile(sww.filename, netcdf_mode_r)
    717 
    718         # Get the variables
    719         x = fid.variables['x']
    720         y = fid.variables['y']
    721         z = fid.variables['elevation']
    722 
    723         volumes = fid.variables['volumes']
    724         time = fid.variables['time']
    725 
    726         # 2D
    727         stage = fid.variables['stage']
    728 
    729         X = x[:]
    730         Y = y[:]
    731         Z = z[:]
    732         V = volumes[:]
    733         T = time[:]
    734         S = stage[:,:]
    735 
    736 #         print "****************************"
    737 #         print "X ",X
    738 #         print "****************************"
    739 #         print "Y ",Y
    740 #         print "****************************"
    741 #         print "Z ",Z
    742 #         print "****************************"
    743 #         print "V ",V
    744 #         print "****************************"
    745 #         print "Time ",T
    746 #         print "****************************"
    747 #         print "Stage ",S
    748 #         print "****************************"
    749 
    750 
    751         fid.close()
    752 
    753         #Cleanup
    754         os.remove(sww.filename)
    755 
    756 
    757 
    758     def test_export_grid(self):
    759         """
    760         test_export_grid(self):
    761         Test that sww information can be converted correctly to asc/prj
    762         format readable by e.g. ArcView
    763         """
    764 
    765         import time, os
    766         from Scientific.IO.NetCDF import NetCDFFile
    767 
    768         try:
    769             os.remove('teg*.sww')
    770         except:
    771             pass
    772 
    773         #Setup
    774         self.domain.set_name('teg')
    775 
    776         prjfile = self.domain.get_name() + '_elevation.prj'
    777         ascfile = self.domain.get_name() + '_elevation.asc'
    778         swwfile = self.domain.get_name() + '.sww'
    779 
    780         self.domain.set_datadir('.')
    781         self.domain.smooth = True
    782         self.domain.set_quantity('elevation', lambda x,y: -x-y)
    783         self.domain.set_quantity('stage', 1.0)
    784 
    785         self.domain.geo_reference = Geo_reference(56,308500,6189000)
    786 
    787         sww = SWW_file(self.domain)
    788         sww.store_connectivity()
    789         sww.store_timestep()
    790         self.domain.evolve_to_end(finaltime = 0.01)
    791         sww.store_timestep()
    792 
    793         cellsize = 0.25
    794         #Check contents
    795         #Get NetCDF
    796 
    797         fid = NetCDFFile(sww.filename, netcdf_mode_r)
    798 
    799         # Get the variables
    800         x = fid.variables['x'][:]
    801         y = fid.variables['y'][:]
    802         z = fid.variables['elevation'][:]
    803         time = fid.variables['time'][:]
    804         stage = fid.variables['stage'][:]
    805 
    806         fid.close()
    807 
    808         #Export to ascii/prj files
    809         sww2dem_batch(self.domain.get_name(),
    810                 quantities = 'elevation',
    811                 cellsize = cellsize,
    812                 verbose = self.verbose,
    813                 format = 'asc')
    814 
    815         #Check asc file
    816         ascid = open(ascfile)
    817         lines = ascid.readlines()
    818         ascid.close()
    819 
    820         L = lines[2].strip().split()
    821         assert L[0].strip().lower() == 'xllcorner'
    822         assert num.allclose(float(L[1].strip().lower()), 308500)
    823 
    824         L = lines[3].strip().split()
    825         assert L[0].strip().lower() == 'yllcorner'
    826         assert num.allclose(float(L[1].strip().lower()), 6189000)
    827 
    828         #Check grid values
    829         for j in range(5):
    830             L = lines[6+j].strip().split()
    831             y = (4-j) * cellsize
    832             for i in range(5):
    833                 assert num.allclose(float(L[i]), -i*cellsize - y)
    834                
    835         #Cleanup
    836         os.remove(prjfile)
    837         os.remove(ascfile)
    838         os.remove(swwfile)
    839 
    840     def test_export_gridII(self):
    841         """
    842         test_export_gridII(self):
    843         Test that sww information can be converted correctly to asc/prj
    844         format readable by e.g. ArcView
    845         """
    846 
    847         import time, os
    848         from Scientific.IO.NetCDF import NetCDFFile
    849 
    850         try:
    851             os.remove('teg*.sww')
    852         except:
    853             pass
    854 
    855         #Setup
    856         self.domain.set_name('tegII')
    857 
    858         swwfile = self.domain.get_name() + '.sww'
    859 
    860         self.domain.set_datadir('.')
    861         self.domain.smooth = True
    862         self.domain.set_quantity('elevation', lambda x,y: -x-y)
    863         self.domain.set_quantity('stage', 1.0)
    864 
    865         self.domain.geo_reference = Geo_reference(56,308500,6189000)
    866 
    867         sww = SWW_file(self.domain)
    868         sww.store_connectivity()
    869         sww.store_timestep()
    870         self.domain.evolve_to_end(finaltime = 0.01)
    871         sww.store_timestep()
    872 
    873         cellsize = 0.25
    874         #Check contents
    875         #Get NetCDF
    876 
    877         fid = NetCDFFile(sww.filename, netcdf_mode_r)
    878 
    879         # Get the variables
    880         x = fid.variables['x'][:]
    881         y = fid.variables['y'][:]
    882         z = fid.variables['elevation'][:]
    883         time = fid.variables['time'][:]
    884         stage = fid.variables['stage'][:]
    885         xmomentum = fid.variables['xmomentum'][:]
    886         ymomentum = fid.variables['ymomentum'][:]       
    887 
    888         #print 'stage', stage
    889         #print 'xmom', xmomentum
    890         #print 'ymom', ymomentum       
    891 
    892         fid.close()
    893 
    894         #Export to ascii/prj files
    895         if True:
    896             sww2dem_batch(self.domain.get_name(),
    897                         quantities = ['elevation', 'depth'],
    898                         cellsize = cellsize,
    899                         verbose = self.verbose,
    900                         format = 'asc')
    901 
    902         else:
    903             sww2dem_batch(self.domain.get_name(),
    904                 quantities = ['depth'],
    905                 cellsize = cellsize,
    906                 verbose = self.verbose,
    907                 format = 'asc')
    908 
    909 
    910             export_grid(self.domain.get_name(),
    911                 quantities = ['elevation'],
    912                 cellsize = cellsize,
    913                 verbose = self.verbose,
    914                 format = 'asc')
    915 
    916         prjfile = self.domain.get_name() + '_elevation.prj'
    917         ascfile = self.domain.get_name() + '_elevation.asc'
    918        
    919         #Check asc file
    920         ascid = open(ascfile)
    921         lines = ascid.readlines()
    922         ascid.close()
    923 
    924         L = lines[2].strip().split()
    925         assert L[0].strip().lower() == 'xllcorner'
    926         assert num.allclose(float(L[1].strip().lower()), 308500)
    927 
    928         L = lines[3].strip().split()
    929         assert L[0].strip().lower() == 'yllcorner'
    930         assert num.allclose(float(L[1].strip().lower()), 6189000)
    931 
    932         #print "ascfile", ascfile
    933         #Check grid values
    934         for j in range(5):
    935             L = lines[6+j].strip().split()
    936             y = (4-j) * cellsize
    937             for i in range(5):
    938                 #print " -i*cellsize - y",  -i*cellsize - y
    939                 #print "float(L[i])", float(L[i])
    940                 assert num.allclose(float(L[i]), -i*cellsize - y)
    941 
    942         #Cleanup
    943         os.remove(prjfile)
    944         os.remove(ascfile)
    945        
    946         #Check asc file
    947         ascfile = self.domain.get_name() + '_depth.asc'
    948         prjfile = self.domain.get_name() + '_depth.prj'
    949         ascid = open(ascfile)
    950         lines = ascid.readlines()
    951         ascid.close()
    952 
    953         L = lines[2].strip().split()
    954         assert L[0].strip().lower() == 'xllcorner'
    955         assert num.allclose(float(L[1].strip().lower()), 308500)
    956 
    957         L = lines[3].strip().split()
    958         assert L[0].strip().lower() == 'yllcorner'
    959         assert num.allclose(float(L[1].strip().lower()), 6189000)
    960 
    961         #Check grid values
    962         for j in range(5):
    963             L = lines[6+j].strip().split()
    964             y = (4-j) * cellsize
    965             for i in range(5):
    966                 #print " -i*cellsize - y",  -i*cellsize - y
    967                 #print "float(L[i])", float(L[i])               
    968                 assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
    969 
    970         #Cleanup
    971         os.remove(prjfile)
    972         os.remove(ascfile)
    973         os.remove(swwfile)
    974 
    975 
    976     def test_export_gridIII(self):
    977         """
    978         test_export_gridIII
    979         Test that sww information can be converted correctly to asc/prj
    980         format readable by e.g. ArcView
    981         """
    982 
    983         import time, os
    984         from Scientific.IO.NetCDF import NetCDFFile
    985 
    986         try:
    987             os.remove('teg*.sww')
    988         except:
    989             pass
    990 
    991         #Setup
    992        
    993         self.domain.set_name('tegIII')
    994 
    995         swwfile = self.domain.get_name() + '.sww'
    996 
    997         self.domain.set_datadir('.')
    998         self.domain.format = 'sww'
    999         self.domain.smooth = True
    1000         self.domain.set_quantity('elevation', lambda x,y: -x-y)
    1001         self.domain.set_quantity('stage', 1.0)
    1002 
    1003         self.domain.geo_reference = Geo_reference(56,308500,6189000)
    1004        
    1005         sww = SWW_file(self.domain)
    1006         sww.store_connectivity()
    1007         sww.store_timestep() #'stage')
    1008         self.domain.evolve_to_end(finaltime = 0.01)
    1009         sww.store_timestep() #'stage')
    1010 
    1011         cellsize = 0.25
    1012         #Check contents
    1013         #Get NetCDF
    1014 
    1015         fid = NetCDFFile(sww.filename, netcdf_mode_r)
    1016 
    1017         # Get the variables
    1018         x = fid.variables['x'][:]
    1019         y = fid.variables['y'][:]
    1020         z = fid.variables['elevation'][:]
    1021         time = fid.variables['time'][:]
    1022         stage = fid.variables['stage'][:]
    1023 
    1024         fid.close()
    1025 
    1026         #Export to ascii/prj files
    1027         extra_name_out = 'yeah'
    1028         if True:
    1029             sww2dem_batch(self.domain.get_name(),
    1030                         quantities = ['elevation', 'depth'],
    1031                         extra_name_out = extra_name_out,
    1032                         cellsize = cellsize,
    1033                         verbose = self.verbose,
    1034                         format = 'asc')
    1035 
    1036         else:
    1037             sww2dem_batch(self.domain.get_name(),
    1038                 quantities = ['depth'],
    1039                 cellsize = cellsize,
    1040                 verbose = self.verbose,
    1041                 format = 'asc')
    1042 
    1043 
    1044             sww2dem_batch(self.domain.get_name(),
    1045                 quantities = ['elevation'],
    1046                 cellsize = cellsize,
    1047                 verbose = self.verbose,
    1048                 format = 'asc')
    1049 
    1050         prjfile = self.domain.get_name() + '_elevation_yeah.prj'
    1051         ascfile = self.domain.get_name() + '_elevation_yeah.asc'
    1052        
    1053         #Check asc file
    1054         ascid = open(ascfile)
    1055         lines = ascid.readlines()
    1056         ascid.close()
    1057 
    1058         L = lines[2].strip().split()
    1059         assert L[0].strip().lower() == 'xllcorner'
    1060         assert num.allclose(float(L[1].strip().lower()), 308500)
    1061 
    1062         L = lines[3].strip().split()
    1063         assert L[0].strip().lower() == 'yllcorner'
    1064         assert num.allclose(float(L[1].strip().lower()), 6189000)
    1065 
    1066         #print "ascfile", ascfile
    1067         #Check grid values
    1068         for j in range(5):
    1069             L = lines[6+j].strip().split()
    1070             y = (4-j) * cellsize
    1071             for i in range(5):
    1072                 #print " -i*cellsize - y",  -i*cellsize - y
    1073                 #print "float(L[i])", float(L[i])
    1074                 assert num.allclose(float(L[i]), -i*cellsize - y)
    1075                
    1076         #Cleanup
    1077         os.remove(prjfile)
    1078         os.remove(ascfile)
    1079        
    1080         #Check asc file
    1081         ascfile = self.domain.get_name() + '_depth_yeah.asc'
    1082         prjfile = self.domain.get_name() + '_depth_yeah.prj'
    1083         ascid = open(ascfile)
    1084         lines = ascid.readlines()
    1085         ascid.close()
    1086 
    1087         L = lines[2].strip().split()
    1088         assert L[0].strip().lower() == 'xllcorner'
    1089         assert num.allclose(float(L[1].strip().lower()), 308500)
    1090 
    1091         L = lines[3].strip().split()
    1092         assert L[0].strip().lower() == 'yllcorner'
    1093         assert num.allclose(float(L[1].strip().lower()), 6189000)
    1094 
    1095         #Check grid values
    1096         for j in range(5):
    1097             L = lines[6+j].strip().split()
    1098             y = (4-j) * cellsize
    1099             for i in range(5):
    1100                 assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
    1101 
    1102         #Cleanup
    1103         os.remove(prjfile)
    1104         os.remove(ascfile)
    1105         os.remove(swwfile)
    1106 
    1107     def test_export_grid_bad(self):
    1108         """Test that sww information can be converted correctly to asc/prj
    1109         format readable by e.g. ArcView
    1110         """
    1111 
    1112         try:
    1113             sww2dem_batch('a_small_round-egg',
    1114                         quantities = ['elevation', 'depth'],
    1115                         cellsize = 99,
    1116                         verbose = self.verbose,
    1117                         format = 'asc')
    1118         except IOError:
    1119             pass
    1120         else:
    1121             self.failUnless(0 ==1,  'Bad input did not throw exception error!')
    1122 
    1123 
    1124641    def test_file_boundary_stsIV_sinewave_ordering(self):
    1125642        """test_file_boundary_stsIV_sinewave_ordering(self):
     
    1195712        xmomentum = fid.variables['xmomentum'][:]
    1196713        ymomentum = fid.variables['ymomentum'][:]
    1197         elevation = fid.variables['elevation'][:]
    1198 
    1199         #print stage
    1200         #print xmomentum
    1201         #print ymomentum
    1202         #print elevation
    1203        
    1204        
     714        elevation = fid.variables['elevation'][:]       
    1205715
    1206716        # Create beginnings of boundary polygon based on sts_boundary
     
    1278788                                                   skip_initial_step=False)):
    1279789            temp_time[i]=domain_time.quantities['stage'].centroid_values[2]
    1280 
    1281 
    1282 
    1283         #print temp_fbound
    1284         #print temp_time
    1285 
    1286         #print domain_fbound.quantities['stage'].vertex_values
    1287         #print domain_time.quantities['stage'].vertex_values
    1288790       
    1289791        assert num.allclose(temp_fbound, temp_time)               
     
    1306808        os.remove(meshname)
    1307809       
    1308        
    1309 
    1310        
    1311        
    1312     def test_file_boundary_sts_time_limit(self):
    1313         """test_file_boundary_stsIV_sinewave_ordering(self):
    1314         Read correct points from ordering file and apply sts to boundary
    1315         This one uses a sine wave and compares to time boundary
    1316        
    1317         This one tests that times used can be limited by upper limit
    1318         """
    1319        
    1320         lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
    1321         bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
    1322         tide = 0.35
    1323         time_step_count = 50
    1324         time_step = 0.1
    1325         times_ref = num.arange(0, time_step_count*time_step, time_step)
    1326        
    1327         n=len(lat_long_points)
    1328         first_tstep=num.ones(n,num.int)
    1329         last_tstep=(time_step_count)*num.ones(n,num.int)
    1330        
    1331         gauge_depth=20*num.ones(n,num.float)
    1332        
    1333         ha1=num.ones((n,time_step_count),num.float)
    1334         ua1=3.*num.ones((n,time_step_count),num.float)
    1335         va1=2.*num.ones((n,time_step_count),num.float)
    1336         for i in range(n):
    1337             ha1[i]=num.sin(times_ref)
    1338        
    1339        
    1340         base_name, files = self.write_mux2(lat_long_points,
    1341                                            time_step_count, time_step,
    1342                                            first_tstep, last_tstep,
    1343                                            depth=gauge_depth,
    1344                                            ha=ha1,
    1345                                            ua=ua1,
    1346                                            va=va1)
    1347 
    1348         # Write order file
    1349         file_handle, order_base_name = tempfile.mkstemp("")
    1350         os.close(file_handle)
    1351         os.remove(order_base_name)
    1352         d=","
    1353         order_file=order_base_name+'order.txt'
    1354         fid=open(order_file,'w')
    1355        
    1356         # Write Header
    1357         header='index, longitude, latitude\n'
    1358         fid.write(header)
    1359         indices=[3,0,1]
    1360         for i in indices:
    1361             line=str(i)+d+str(lat_long_points[i][1])+d+\
    1362                 str(lat_long_points[i][0])+"\n"
    1363             fid.write(line)
    1364         fid.close()
    1365 
    1366         sts_file=base_name
    1367         urs2sts(base_name, basename_out=sts_file,
    1368                 ordering_filename=order_file,
    1369                 mean_stage=tide,
    1370                 verbose=False)
    1371         self.delete_mux(files)
    1372        
    1373        
    1374        
    1375         # Now read the sts file and check that values have been stored correctly.
    1376         fid = NetCDFFile(sts_file + '.sts')
    1377 
    1378         # Check the time vector
    1379         times = fid.variables['time'][:]
    1380         starttime = fid.starttime[0]
    1381         #print times
    1382         #print starttime
    1383 
    1384         # Check sts quantities
    1385         stage = fid.variables['stage'][:]
    1386         xmomentum = fid.variables['xmomentum'][:]
    1387         ymomentum = fid.variables['ymomentum'][:]
    1388         elevation = fid.variables['elevation'][:]
    1389 
    1390        
    1391 
    1392         # Create beginnings of boundary polygon based on sts_boundary
    1393         boundary_polygon = create_sts_boundary(base_name)
    1394        
    1395         os.remove(order_file)
    1396 
    1397         # Append the remaining part of the boundary polygon to be defined by
    1398         # the user
    1399         bounding_polygon_utm=[]
    1400         for point in bounding_polygon:
    1401             zone,easting,northing=redfearn(point[0],point[1])
    1402             bounding_polygon_utm.append([easting,northing])
    1403 
    1404         boundary_polygon.append(bounding_polygon_utm[3])
    1405         boundary_polygon.append(bounding_polygon_utm[4])
    1406 
    1407         #print 'boundary_polygon', boundary_polygon
    1408        
    1409 
    1410         assert num.allclose(bounding_polygon_utm,boundary_polygon)
    1411 
    1412 
    1413         extent_res=1000000
    1414         meshname = 'urs_test_mesh' + '.tsh'
    1415         interior_regions=None
    1416         boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
    1417        
    1418         # have to change boundary tags from last example because now bounding
    1419         # polygon starts in different place.
    1420         create_mesh_from_regions(boundary_polygon,
    1421                                  boundary_tags=boundary_tags,
    1422                                  maximum_triangle_area=extent_res,
    1423                                  filename=meshname,
    1424                                  interior_regions=interior_regions,
    1425                                  verbose=False)
    1426        
    1427         domain_fbound = Domain(meshname)
    1428         domain_fbound.set_quantity('stage', tide)
    1429        
    1430        
    1431         Bf = File_boundary(sts_file+'.sts',
    1432                            domain_fbound,
    1433                            boundary_polygon=boundary_polygon)
    1434         time_vec = Bf.F.get_time()
    1435         assert num.allclose(Bf.F.starttime, starttime)
    1436         assert num.allclose(time_vec, times_ref)                                   
    1437        
    1438         for time_limit in [0.1, 0.2, 0.5, 1.0, 2.2, 3.0, 4.3, 6.0, 10.0]:
    1439             Bf = File_boundary(sts_file+'.sts',
    1440                                domain_fbound,
    1441                                time_limit=time_limit+starttime,
    1442                                boundary_polygon=boundary_polygon)
    1443        
    1444             time_vec = Bf.F.get_time()
    1445             assert num.allclose(Bf.F.starttime, starttime)           
    1446             assert num.alltrue(time_vec < time_limit)
    1447            
    1448            
    1449         try:   
    1450             Bf = File_boundary(sts_file+'.sts',
    1451                                domain_fbound,
    1452                                time_limit=-1+starttime,
    1453                                boundary_polygon=boundary_polygon)           
    1454             time_vec = Bf.F.get_time()   
    1455             print time_vec   
    1456         except AssertionError:
    1457             pass
    1458         else:
    1459             raise Exception, 'Should have raised Exception here'
    1460 
    1461     #### END TESTS FOR URS 2 SWW  ###
    1462 
    1463 
    1464     def test_triangulation(self):
    1465         #
    1466         # 
    1467        
    1468         filename = tempfile.mktemp("_data_manager.sww")
    1469         outfile = NetCDFFile(filename, netcdf_mode_w)
    1470         points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
    1471         volumes = (0,1,2)
    1472         elevation = [0,1,2]
    1473         new_origin = None
    1474         new_origin = Geo_reference(56, 0, 0)
    1475         times = [0, 10]
    1476         number_of_volumes = len(volumes)
    1477         number_of_points = len(points_utm)
    1478         sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
    1479         sww.store_header(outfile, times, number_of_volumes,
    1480                          number_of_points, description='fully sick testing',
    1481                          verbose=self.verbose,sww_precision=netcdf_float)
    1482         sww.store_triangulation(outfile, points_utm, volumes,
    1483                                 elevation,  new_origin=new_origin,
    1484                                 verbose=self.verbose)       
    1485         outfile.close()
    1486         fid = NetCDFFile(filename)
    1487 
    1488         x = fid.variables['x'][:]
    1489         y = fid.variables['y'][:]
    1490         fid.close()
    1491 
    1492         assert num.allclose(num.array(map(None, x,y)), points_utm)
    1493         os.remove(filename)
    1494 
    1495        
    1496     def test_triangulationII(self):
    1497         #
    1498         # 
    1499        
    1500         filename = tempfile.mktemp("_data_manager.sww")
    1501         outfile = NetCDFFile(filename, netcdf_mode_w)
    1502         points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
    1503         volumes = (0,1,2)
    1504         elevation = [0,1,2]
    1505         new_origin = None
    1506         #new_origin = Geo_reference(56, 0, 0)
    1507         times = [0, 10]
    1508         number_of_volumes = len(volumes)
    1509         number_of_points = len(points_utm)
    1510         sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
    1511         sww.store_header(outfile, times, number_of_volumes,
    1512                          number_of_points, description='fully sick testing',
    1513                          verbose=self.verbose,sww_precision=netcdf_float)
    1514         sww.store_triangulation(outfile, points_utm, volumes,
    1515                                 new_origin=new_origin,
    1516                                 verbose=self.verbose)
    1517         sww.store_static_quantities(outfile, elevation=elevation)                               
    1518                                
    1519         outfile.close()
    1520         fid = NetCDFFile(filename)
    1521 
    1522         x = fid.variables['x'][:]
    1523         y = fid.variables['y'][:]
    1524         results_georef = Geo_reference()
    1525         results_georef.read_NetCDF(fid)
    1526         assert results_georef == Geo_reference(DEFAULT_ZONE, 0, 0)
    1527         fid.close()
    1528 
    1529         assert num.allclose(num.array(map(None, x,y)), points_utm)
    1530         os.remove(filename)
    1531 
    1532        
    1533     def test_triangulation_new_origin(self):
    1534         #
    1535         # 
    1536        
    1537         filename = tempfile.mktemp("_data_manager.sww")
    1538         outfile = NetCDFFile(filename, netcdf_mode_w)
    1539         points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
    1540         volumes = (0,1,2)
    1541         elevation = [0,1,2]
    1542         new_origin = None
    1543         new_origin = Geo_reference(56, 1, 554354)
    1544         points_utm = new_origin.change_points_geo_ref(points_utm)
    1545         times = [0, 10]
    1546         number_of_volumes = len(volumes)
    1547         number_of_points = len(points_utm)
    1548         sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
    1549         sww.store_header(outfile, times, number_of_volumes,
    1550                          number_of_points, description='fully sick testing',
    1551                          verbose=self.verbose,sww_precision=netcdf_float)
    1552         sww.store_triangulation(outfile, points_utm, volumes,
    1553                                 elevation,  new_origin=new_origin,
    1554                                 verbose=self.verbose)
    1555         outfile.close()
    1556         fid = NetCDFFile(filename)
    1557 
    1558         x = fid.variables['x'][:]
    1559         y = fid.variables['y'][:]
    1560         results_georef = Geo_reference()
    1561         results_georef.read_NetCDF(fid)
    1562         assert results_georef == new_origin
    1563         fid.close()
    1564 
    1565         absolute = Geo_reference(56, 0,0)
    1566         assert num.allclose(num.array(
    1567             absolute.change_points_geo_ref(map(None, x,y),
    1568                                            new_origin)),points_utm)
    1569        
    1570         os.remove(filename)
    1571        
    1572     def test_triangulation_points_georeference(self):
    1573         #
    1574         # 
    1575        
    1576         filename = tempfile.mktemp("_data_manager.sww")
    1577         outfile = NetCDFFile(filename, netcdf_mode_w)
    1578         points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
    1579         volumes = (0,1,2)
    1580         elevation = [0,1,2]
    1581         new_origin = None
    1582         points_georeference = Geo_reference(56, 1, 554354)
    1583         points_utm = points_georeference.change_points_geo_ref(points_utm)
    1584         times = [0, 10]
    1585         number_of_volumes = len(volumes)
    1586         number_of_points = len(points_utm)
    1587         sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
    1588         sww.store_header(outfile, times, number_of_volumes,
    1589                          number_of_points, description='fully sick testing',
    1590                          verbose=self.verbose,sww_precision=netcdf_float)
    1591         sww.store_triangulation(outfile, points_utm, volumes,
    1592                                 elevation,  new_origin=new_origin,
    1593                                 points_georeference=points_georeference,
    1594                                 verbose=self.verbose)       
    1595         outfile.close()
    1596         fid = NetCDFFile(filename)
    1597 
    1598         x = fid.variables['x'][:]
    1599         y = fid.variables['y'][:]
    1600         results_georef = Geo_reference()
    1601         results_georef.read_NetCDF(fid)
    1602         assert results_georef == points_georeference
    1603         fid.close()
    1604 
    1605         assert num.allclose(num.array(map(None, x,y)), points_utm)
    1606         os.remove(filename)
    1607        
    1608     def test_triangulation_2_geo_refs(self):
    1609         #
    1610         # 
    1611        
    1612         filename = tempfile.mktemp("_data_manager.sww")
    1613         outfile = NetCDFFile(filename, netcdf_mode_w)
    1614         points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
    1615         volumes = (0,1,2)
    1616         elevation = [0,1,2]
    1617         new_origin = Geo_reference(56, 1, 1)
    1618         points_georeference = Geo_reference(56, 0, 0)
    1619         points_utm = points_georeference.change_points_geo_ref(points_utm)
    1620         times = [0, 10]
    1621         number_of_volumes = len(volumes)
    1622         number_of_points = len(points_utm)
    1623         sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])       
    1624         sww.store_header(outfile, times, number_of_volumes,
    1625                          number_of_points, description='fully sick testing',
    1626                          verbose=self.verbose,sww_precision=netcdf_float)
    1627         sww.store_triangulation(outfile, points_utm, volumes,
    1628                                 elevation,  new_origin=new_origin,
    1629                                 points_georeference=points_georeference,
    1630                                 verbose=self.verbose)       
    1631         outfile.close()
    1632         fid = NetCDFFile(filename)
    1633 
    1634         x = fid.variables['x'][:]
    1635         y = fid.variables['y'][:]
    1636         results_georef = Geo_reference()
    1637         results_georef.read_NetCDF(fid)
    1638         assert results_georef == new_origin
    1639         fid.close()
    1640 
    1641 
    1642         absolute = Geo_reference(56, 0,0)
    1643         assert num.allclose(num.array(
    1644             absolute.change_points_geo_ref(map(None, x,y),
    1645                                            new_origin)),points_utm)
    1646         os.remove(filename)
    1647 
    1648810 
    1649811    def test_points2polygon(self): 
     
    1664826        assert (polygon == [[0.0, 0.0],[1.0, 0.0],[0.0, 1.0]])
    1665827       
    1666    
    1667     def test_csv2polygons(self):
    1668         """test_csv2polygons
    1669         """
    1670        
    1671         path = get_pathname_from_package('anuga.shallow_water')               
    1672         testfile = os.path.join(path, 'polygon_values_example.csv')               
    1673 
    1674         polygons, values = load_csv_as_polygons(testfile,
    1675                                         value_name='floors')
    1676 
    1677         assert len(polygons) == 7, 'Must have seven polygons'
    1678         assert len(values) == 7, 'Must have seven values'
    1679            
    1680         # Known floor values
    1681         floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
    1682        
    1683         # Known polygon values
    1684         known_polys = {'1': [[422681.61,871117.55],
    1685                              [422691.02,871117.60],
    1686                              [422690.87,871084.23],
    1687                              [422649.36,871081.85],
    1688                              [422649.36,871080.39],
    1689                              [422631.86,871079.50],
    1690                              [422631.72,871086.75],
    1691                              [422636.75,871087.20],
    1692                              [422636.75,871091.50],
    1693                              [422649.66,871092.09],
    1694                              [422649.83,871084.91],
    1695                              [422652.94,871084.90],
    1696                              [422652.84,871092.39],
    1697                              [422681.83,871093.73],
    1698                              [422681.61,871117.55]],
    1699                        '2': [[422664.22,870785.46],
    1700                              [422672.48,870780.14],
    1701                              [422668.17,870772.62],
    1702                              [422660.35,870777.17],
    1703                              [422664.22,870785.46]],
    1704                        '3': [[422661.30,871215.06],
    1705                              [422667.50,871215.70],
    1706                              [422668.30,871204.86],
    1707                              [422662.21,871204.33],
    1708                              [422661.30,871215.06]],
    1709                        '4': [[422473.44,871191.22],
    1710                              [422478.33,871192.26],
    1711                              [422479.52,871186.03],
    1712                              [422474.78,871185.14],
    1713                              [422473.44,871191.22]],
    1714                        '5': [[422369.69,871049.29],
    1715                              [422378.63,871053.58],
    1716                              [422383.91,871044.51],
    1717                              [422374.97,871040.32],
    1718                              [422369.69,871049.29]],
    1719                        '8': [[422730.56,871203.13],
    1720                              [422734.10,871204.90],
    1721                              [422735.26,871202.18],
    1722                              [422731.87,871200.58],
    1723                              [422730.56,871203.13]],
    1724                        '9': [[422659.85,871213.80],
    1725                              [422660.91,871210.97],
    1726                              [422655.42,871208.85],
    1727                              [422654.36,871211.68],
    1728                              [422659.85,871213.80]]
    1729                        }       
    1730        
    1731 
    1732        
    1733                
    1734         for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
    1735             assert id in polygons.keys()
    1736             assert id in values.keys()           
    1737 
    1738             assert int(values[id]) == int(floors[id])
    1739             assert len(polygons[id]) == len(known_polys[id])
    1740             assert num.allclose(polygons[id], known_polys[id])
    1741 
    1742 
    1743     def test_csv2polygons_with_clipping(self):
    1744         """test_csv2polygons with optional clipping
    1745         """
    1746         #FIXME(Ole): Not Done!!
    1747        
    1748         path = get_pathname_from_package('anuga.shallow_water')               
    1749         testfile = os.path.join(path, 'polygon_values_example.csv')               
    1750 
    1751         polygons, values = load_csv_as_polygons(testfile,
    1752                                         value_name='floors',
    1753                                         clipping_polygons=None)
    1754 
    1755         assert len(polygons) == 7, 'Must have seven polygons'
    1756         assert len(values) == 7, 'Must have seven values'
    1757            
    1758         # Known floor values
    1759         floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
    1760        
    1761         # Known polygon values
    1762         known_polys = {'1': [[422681.61,871117.55],
    1763                              [422691.02,871117.60],
    1764                              [422690.87,871084.23],
    1765                              [422649.36,871081.85],
    1766                              [422649.36,871080.39],
    1767                              [422631.86,871079.50],
    1768                              [422631.72,871086.75],
    1769                              [422636.75,871087.20],
    1770                              [422636.75,871091.50],
    1771                              [422649.66,871092.09],
    1772                              [422649.83,871084.91],
    1773                              [422652.94,871084.90],
    1774                              [422652.84,871092.39],
    1775                              [422681.83,871093.73],
    1776                              [422681.61,871117.55]],
    1777                        '2': [[422664.22,870785.46],
    1778                              [422672.48,870780.14],
    1779                              [422668.17,870772.62],
    1780                              [422660.35,870777.17],
    1781                              [422664.22,870785.46]],
    1782                        '3': [[422661.30,871215.06],
    1783                              [422667.50,871215.70],
    1784                              [422668.30,871204.86],
    1785                              [422662.21,871204.33],
    1786                              [422661.30,871215.06]],
    1787                        '4': [[422473.44,871191.22],
    1788                              [422478.33,871192.26],
    1789                              [422479.52,871186.03],
    1790                              [422474.78,871185.14],
    1791                              [422473.44,871191.22]],
    1792                        '5': [[422369.69,871049.29],
    1793                              [422378.63,871053.58],
    1794                              [422383.91,871044.51],
    1795                              [422374.97,871040.32],
    1796                              [422369.69,871049.29]],
    1797                        '8': [[422730.56,871203.13],
    1798                              [422734.10,871204.90],
    1799                              [422735.26,871202.18],
    1800                              [422731.87,871200.58],
    1801                              [422730.56,871203.13]],
    1802                        '9': [[422659.85,871213.80],
    1803                              [422660.91,871210.97],
    1804                              [422655.42,871208.85],
    1805                              [422654.36,871211.68],
    1806                              [422659.85,871213.80]]
    1807                        }       
    1808        
    1809 
    1810        
    1811                
    1812         for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
    1813             assert id in polygons.keys()
    1814             assert id in values.keys()           
    1815 
    1816             assert int(values[id]) == int(floors[id])
    1817             assert len(polygons[id]) == len(known_polys[id])
    1818             assert num.allclose(polygons[id], known_polys[id])
    1819 
    1820 
    1821 
    1822 
    1823    
    1824     def test_csv2building_polygons(self):
    1825         """test_csv2building_polygons
    1826         """
    1827        
    1828         path = get_pathname_from_package('anuga.shallow_water')               
    1829         testfile = os.path.join(path, 'polygon_values_example.csv')               
    1830 
    1831         polygons, values = load_csv_as_building_polygons(testfile,
    1832                                                  floor_height=3)
    1833 
    1834         assert len(polygons) == 7, 'Must have seven polygons'
    1835         assert len(values) == 7, 'Must have seven values'
    1836            
    1837         # Known floor values
    1838         floors = {'1': 6, '2': 0, '3': 3, '4': 6, '5': 0, '8': 3, '9': 3}
    1839        
    1840                
    1841         for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
    1842             assert id in polygons.keys()
    1843             assert id in values.keys()           
    1844 
    1845             assert float(values[id]) == float(floors[id])
    1846828
    1847829
  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_basic.py

    r7575 r7866  
    66import tempfile
    77
     8import anuga
     9
    810from anuga.config import g, epsilon
    911from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    1012from anuga.utilities.numerical_tools import mean
    11 from anuga.utilities.polygon import is_inside_polygon
     13from anuga.geometry.polygon import is_inside_polygon
    1214from anuga.coordinate_transforms.geo_reference import Geo_reference
    1315from anuga.abstract_2d_finite_volumes.quantity import Quantity
     
    20852087        can be controlled this way.
    20862088        """
    2087 
    2088         #---------------------------------------------------------------------
    2089         # Import necessary modules
    2090         #---------------------------------------------------------------------
    2091         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    2092         from anuga.shallow_water import Domain
    2093         from anuga.shallow_water import Reflective_boundary
    2094         from anuga.shallow_water import Dirichlet_boundary
    2095         from anuga.shallow_water import Time_boundary
    2096 
     2089       
    20972090        #---------------------------------------------------------------------
    20982091        # Setup computational domain
     
    21082101                                                       len1=length,
    21092102                                                       len2=width)
    2110         domain = Domain(points, vertices, boundary)
     2103        domain = anuga.Domain(points, vertices, boundary)
    21112104        domain.set_name('channel_variable_test')  # Output name
    21122105        domain.set_quantities_to_be_stored({'elevation': 2,
  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py

    r7616 r7866  
    99from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    1010from anuga.utilities.numerical_tools import mean
    11 from anuga.utilities.polygon import is_inside_polygon
     11from anuga.geometry.polygon import is_inside_polygon
    1212from anuga.coordinate_transforms.geo_reference import Geo_reference
    1313from anuga.abstract_2d_finite_volumes.quantity import Quantity
  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_conservation.py

    r7733 r7866  
    66import tempfile
    77
     8import anuga
     9
    810from anuga.config import g, epsilon
    911from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    1012from anuga.utilities.numerical_tools import mean
    11 from anuga.utilities.polygon import is_inside_polygon
     13from anuga.geometry.polygon import is_inside_polygon
    1214from anuga.coordinate_transforms.geo_reference import Geo_reference
    1315from anuga.abstract_2d_finite_volumes.quantity import Quantity
     
    6062
    6163        # Boundary conditions (reflective everywhere)
    62         Br = Reflective_boundary(domain)
     64        Br = anuga.Reflective_boundary(domain)
    6365        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    6466
     
    8890        """
    8991
    90         from mesh_factory import rectangular
    91 
    9292        # Create basic mesh
    93         points, vertices, boundary = rectangular(6, 6)
     93        points, vertices, boundary = anuga.rectangular(6, 6)
    9494
    9595        # Create shallow water domain
    96         domain = Domain(points, vertices, boundary)
     96        domain = anuga.Domain(points, vertices, boundary)
    9797        domain.smooth = False
    9898        domain.default_order = 2
     
    107107
    108108        # Boundary conditions (reflective everywhere)
    109         Br = Reflective_boundary(domain)
     109        Br = anuga.Reflective_boundary(domain)
    110110        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    111111
     
    140140
    141141        # Create shallow water domain
    142         domain = Domain(points, vertices, boundary)
     142        domain = anuga.Domain(points, vertices, boundary)
    143143        domain.smooth = False
    144144
     
    163163
    164164        # Boundary conditions (reflective everywhere)
    165         Br = Reflective_boundary(domain)
     165        Br = anuga.Reflective_boundary(domain)
    166166        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    167167
     
    200200        """
    201201
    202         from mesh_factory import rectangular
    203 
    204202        # Create basic mesh
    205         points, vertices, boundary = rectangular(6, 6)
     203        points, vertices, boundary = anuga.rectangular(6, 6)
    206204
    207205        # Create shallow water domain
    208         domain = Domain(points, vertices, boundary)
     206        domain = anuga.Domain(points, vertices, boundary)
    209207        domain.smooth = False
    210208        domain.default_order = 2
     
    232230
    233231        # Boundary conditions (reflective everywhere)
    234         Br = Reflective_boundary(domain)
     232        Br = anuga.Reflective_boundary(domain)
    235233        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    236234
     
    271269        """
    272270
    273         from mesh_factory import rectangular
    274 
    275271        # Create basic mesh
    276         points, vertices, boundary = rectangular(6, 6)
     272        points, vertices, boundary = anuga.rectangular(6, 6)
    277273
    278274        # Create shallow water domain
    279         domain = Domain(points, vertices, boundary)
     275        domain = anuga.Domain(points, vertices, boundary)
    280276        domain.smooth = False
    281277        domain.default_order = 2
     
    290286
    291287        # Boundary conditions (reflective everywhere)
    292         Br = Reflective_boundary(domain)
    293         Bleft = Dirichlet_boundary([0.5, 0, 0])
    294         Bright = Dirichlet_boundary([0.1, 0, 0])
     288        Br = anuga.Reflective_boundary(domain)
     289        Bleft = anuga.Dirichlet_boundary([0.5, 0, 0])
     290        Bright = anuga.Dirichlet_boundary([0.1, 0, 0])
    295291        domain.set_boundary({'left': Bleft, 'right': Bright,
    296292                             'top': Br, 'bottom': Br})
     
    362358
    363359        # Boundaries
    364         R = Reflective_boundary(domain)
     360        R = anuga.Reflective_boundary(domain)
    365361        domain.set_boundary({'left': R, 'right': R, 'top':R, 'bottom': R})
    366362
     
    384380        Test that total volume can be computed correctly
    385381        """           
    386 
    387         #----------------------------------------------------------------------
    388         # Import necessary modules
    389         #----------------------------------------------------------------------
    390         from anuga.abstract_2d_finite_volumes.mesh_factory \
    391                 import rectangular_cross
    392         from anuga.shallow_water import Domain
    393382
    394383        #----------------------------------------------------------------------
     
    449438        verbose = False
    450439
    451         #----------------------------------------------------------------------
    452         # Import necessary modules
    453         #----------------------------------------------------------------------
    454 
    455         from anuga.abstract_2d_finite_volumes.mesh_factory \
    456                 import rectangular_cross
    457         from anuga.shallow_water import Domain
    458         from anuga.shallow_water.shallow_water_domain import Reflective_boundary
    459         from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
    460         from anuga.shallow_water.forcing import Inflow
    461         from anuga.shallow_water.data_manager \
    462                 import get_flow_through_cross_section
     440
    463441
    464442        #----------------------------------------------------------------------
     
    500478        #----------------------------------------------------------------------
    501479
    502         Br = Reflective_boundary(domain)      # Solid reflective wall
     480        Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
    503481               
    504482        # Constant flow in and out of domain
    505483        # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s
    506         Bi = Dirichlet_boundary([d, uh, vh])
    507         Bo = Dirichlet_boundary([d, uh, vh])
     484        Bi = anuga.Dirichlet_boundary([d, uh, vh])
     485        Bo = anuga.Dirichlet_boundary([d, uh, vh])
    508486
    509487        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
     
    544522       
    545523
    546         #---------------------------------------------------------------------
    547         # Import necessary modules
    548         #---------------------------------------------------------------------
    549         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    550         from anuga.shallow_water import Domain
    551         from anuga.shallow_water.shallow_water_domain import Reflective_boundary
    552         from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
    553         from anuga.shallow_water.forcing import Inflow
    554         from anuga.shallow_water.data_manager import get_flow_through_cross_section
    555 
    556524        #----------------------------------------------------------------------
    557525        # Setup computational domain
     
    569537
    570538
    571         domain = Domain(points, vertices, boundary)   
     539        domain = anuga.Domain(points, vertices, boundary)   
    572540        domain.set_name('Inflow_volume_test')              # Output name
    573541               
     
    593561
    594562        # Fixed Flowrate onto Area
    595         fixed_inflow = Inflow(domain,
     563        fixed_inflow = anuga.Inflow(domain,
    596564                              center=(10.0, 10.0),
    597565                              radius=5.00,
     
    604572        #----------------------------------------------------------------------
    605573
    606         Br = Reflective_boundary(domain) # Solid reflective wall
     574        Br = anuga.Reflective_boundary(domain) # Solid reflective wall
    607575        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    608576
     
    646614       
    647615
    648         #---------------------------------------------------------------------
    649         # Import necessary modules
    650         #---------------------------------------------------------------------
    651         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    652         from anuga.shallow_water import Domain
    653         from anuga.shallow_water.shallow_water_domain import Reflective_boundary
    654         from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
    655         from anuga.shallow_water.shallow_water_domain import Rainfall
    656         from anuga.shallow_water.data_manager import get_flow_through_cross_section
    657 
    658616        #----------------------------------------------------------------------
    659617        # Setup computational domain
     
    666624       
    667625
    668         points, vertices, boundary = rectangular_cross(int(length/dx),
     626        points, vertices, boundary = anuga.rectangular_cross(int(length/dx),
    669627                                                       int(width/dy),
    670628                                                       len1=length, len2=width)
    671629
    672630
    673         domain = Domain(points, vertices, boundary)   
     631        domain = anuga.Domain(points, vertices, boundary)   
    674632        domain.set_name('Rain_volume_test')              # Output name
    675633               
     
    695653
    696654        # Fixed rain onto small circular area
    697         fixed_rain = Rainfall(domain,
     655        fixed_rain = anuga.Rainfall(domain,
    698656                              center=(10.0, 10.0),
    699657                              radius=5.00,
     
    706664        #----------------------------------------------------------------------
    707665
    708         Br = Reflective_boundary(domain) # Solid reflective wall
     666        Br = anuga.Reflective_boundary(domain) # Solid reflective wall
    709667        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    710668
     
    755713       
    756714
    757         #---------------------------------------------------------------------
    758         # Import necessary modules
    759         #---------------------------------------------------------------------
    760         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    761         from anuga.shallow_water import Domain
    762         from anuga.shallow_water.shallow_water_domain import Reflective_boundary
    763         from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
    764         from anuga.shallow_water.shallow_water_domain import Rainfall
    765         from anuga.shallow_water.data_manager import get_flow_through_cross_section
    766 
    767715        #----------------------------------------------------------------------
    768716        # Setup computational domain
     
    815763        #----------------------------------------------------------------------
    816764
    817         Br = Reflective_boundary(domain) # Solid reflective wall
    818         Bt = Transmissive_stage_zero_momentum_boundary(domain)
    819         Bd = Dirichlet_boundary([-10, 0, 0])
     765        Br = anuga.Reflective_boundary(domain) # Solid reflective wall
     766        Bt = anuga.Transmissive_stage_zero_momentum_boundary(domain)
     767        Bd = anuga.Dirichlet_boundary([-10, 0, 0])
    820768        domain.set_boundary({'left': Bt, 'right': Bd, 'top': Bt, 'bottom': Bt})
    821769
  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_cross_sections.py

    r7559 r7866  
    66import tempfile
    77
     8import anuga
     9
    810from anuga.config import g, epsilon
    911from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    1012from anuga.utilities.numerical_tools import mean
    11 from anuga.utilities.polygon import is_inside_polygon
     13from anuga.geometry.polygon import is_inside_polygon
    1214from anuga.coordinate_transforms.geo_reference import Geo_reference
    13 from anuga.abstract_2d_finite_volumes.quantity import Quantity
    14 from anuga.geospatial_data.geospatial_data import Geospatial_data
    15 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    1615
    1716from anuga.utilities.system_tools import get_pathname_from_package
     
    7675        uh = u*h
    7776
    78         Br = Reflective_boundary(domain)     # Side walls
    79         Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
     77        Br = anuga.Reflective_boundary(domain)     # Side walls
     78        Bd = anuga.Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
    8079
    8180
     
    161160        uh = u*h
    162161
    163         Br = Reflective_boundary(domain)       # Side walls
    164         Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
     162        Br = anuga.Reflective_boundary(domain)     # Side walls
     163        Bd = anuga.Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
    165164
    166165        # Initial conditions
     
    252251        uh = u*h
    253252
    254         Br = Reflective_boundary(domain)       # Side walls
    255         Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
     253        Br = anuga.Reflective_boundary(domain)     # Side walls
     254        Bd = anuga.Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
    256255
    257256        # Initial conditions
  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_distribute.py

    r7573 r7866  
    99from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    1010from anuga.utilities.numerical_tools import mean
    11 from anuga.utilities.polygon import is_inside_polygon
     11from anuga.geometry.polygon import is_inside_polygon
    1212from anuga.coordinate_transforms.geo_reference import Geo_reference
    1313from anuga.abstract_2d_finite_volumes.quantity import Quantity
  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_forcing_terms.py

    r7733 r7866  
    55from math import pi, sqrt
    66import tempfile
     7import anuga
    78
    89from anuga.config import g, epsilon
    910from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    1011from anuga.utilities.numerical_tools import mean
    11 from anuga.utilities.polygon import is_inside_polygon
     12from anuga.geometry.polygon import is_inside_polygon
    1213from anuga.coordinate_transforms.geo_reference import Geo_reference
    1314from anuga.abstract_2d_finite_volumes.quantity import Quantity
     
    394395        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    395396
    396         domain = Domain(points, vertices)
     397        domain = anuga.Domain(points, vertices)
    397398
    398399        #Flat surface with 1m of water
     
    401402        domain.set_quantity('friction', 0)
    402403
    403         Br = Reflective_boundary(domain)
     404        Br = anuga.Reflective_boundary(domain)
    404405        domain.set_boundary({'exterior': Br})
    405406
     
    408409        phi = 135
    409410        domain.forcing_terms = []
    410         domain.forcing_terms.append(Wind_stress(s, phi))
     411        domain.forcing_terms.append(anuga.Wind_stress(s, phi))
    411412
    412413        domain.compute_forcing_terms()
     
    450451        domain.set_quantity('friction', 0)
    451452
    452         Br = Reflective_boundary(domain)
     453        Br = anuga.Reflective_boundary(domain)
    453454        domain.set_boundary({'exterior': Br})
    454455
     
    459460        phi = 135
    460461        domain.forcing_terms = []
    461         domain.forcing_terms.append(Wind_stress(s=speed, phi=angle))
     462        domain.forcing_terms.append(anuga.Wind_stress(s=speed, phi=angle))
    462463
    463464        domain.compute_forcing_terms()
     
    522523        domain.set_quantity('friction', 0)
    523524
    524         Br = Reflective_boundary(domain)
     525        Br = anuga.Reflective_boundary(domain)
    525526        domain.set_boundary({'exterior': Br})
    526527
     
    543544        fid.close()
    544545
    545         # Convert ASCII file to NetCDF (Which is what we really like!)
    546         from data_manager import timefile2netcdf
    547 
    548         timefile2netcdf(filename)
     546        anuga.timefile2netcdf(filename+'.txt', filename+'.tms')
    549547        os.remove(filename + '.txt')
    550548
     
    616614        domain.set_quantity('friction', 0)
    617615
    618         Br = Reflective_boundary(domain)
     616        Br = anuga.Reflective_boundary(domain)
    619617        domain.set_boundary({'exterior': Br})
    620618
     
    623621        # Write wind stress file (ensure that domain.time is covered)
    624622        # Take x=1 and y=0
    625         filename = 'test_windstress_from_file'
     623        filename = 'test_windstress_from_file.txt'
     624        file_out = 'test_windstress_from_file.tms'
    626625        start = time.mktime(time.strptime('2000', '%Y'))
    627         fid = open(filename + '.txt', 'w')
     626        fid = open(filename, 'w')
    628627        dt = 0.5    # Half second interval
    629628        t = 0.0
     
    635634        fid.close()
    636635
    637         # Convert ASCII file to NetCDF (Which is what we really like!)
    638         from data_manager import timefile2netcdf
    639 
    640         timefile2netcdf(filename, time_as_seconds=True)
    641         os.remove(filename + '.txt')
     636        anuga.timefile2netcdf(filename, file_out, time_as_seconds=True)
     637        os.remove(filename)
    642638
    643639        # Setup wind stress
    644         F = file_function(filename + '.tms',
     640        F = file_function(file_out,
    645641                          quantities=['Attribute0', 'Attribute1'])
    646         os.remove(filename + '.tms')
    647 
    648         W = Wind_stress(F)
     642        os.remove(file_out)
     643
     644        W = anuga.Wind_stress(F)
    649645
    650646        domain.forcing_terms = []
     
    709705        domain.set_quantity('friction', 0)
    710706
    711         Br = Reflective_boundary(domain)
     707        Br = anuga.Reflective_boundary(domain)
    712708        domain.set_boundary({'exterior': Br})
    713709
     
    718714
    719715        try:
    720             domain.forcing_terms.append(Wind_stress(s=scalar_func_list,
     716            domain.forcing_terms.append(anuga.Wind_stress(s=scalar_func_list,
    721717                                                    phi=angle))
    722718        except AssertionError:
     
    763759        domain.set_quantity('friction', 0)
    764760
    765         Br = Reflective_boundary(domain)
     761        Br = anuga.Reflective_boundary(domain)
    766762        domain.set_boundary({'exterior': Br})
    767763
    768764        # Setup only one forcing term, constant rainfall
    769765        domain.forcing_terms = []
    770         domain.forcing_terms.append(Rainfall(domain, rate=2.0))
     766        domain.forcing_terms.append(anuga.Rainfall(domain, rate=2.0))
    771767
    772768        domain.compute_forcing_terms()
     
    795791        domain.set_quantity('friction', 0)
    796792
    797         Br = Reflective_boundary(domain)
     793        Br = anuga.Reflective_boundary(domain)
    798794        domain.set_boundary({'exterior': Br})
    799795
     
    801797        # restricted to a polygon enclosing triangle #1 (bce)
    802798        domain.forcing_terms = []
    803         R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
     799        R = anuga.Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
    804800
    805801        assert num.allclose(R.exchange_area, 2)
     
    826822        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    827823
    828         domain = Domain(points, vertices)
     824        domain = anuga.Domain(points, vertices)
    829825
    830826        # Flat surface with 1m of water
     
    833829        domain.set_quantity('friction', 0)
    834830
    835         Br = Reflective_boundary(domain)
     831        Br = anuga.Reflective_boundary(domain)
    836832        domain.set_boundary({'exterior': Br})
    837833
     
    839835        # restricted to a polygon enclosing triangle #1 (bce)
    840836        domain.forcing_terms = []
    841         R = Rainfall(domain,
     837        R = anuga.Rainfall(domain,
    842838                     rate=lambda t: 3*t + 7,
    843839                     polygon = [[1,1], [2,1], [2,2], [1,2]])
     
    870866        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    871867
    872         domain = Domain(points, vertices)
     868        domain = anuga.Domain(points, vertices)
    873869
    874870        # Flat surface with 1m of water
     
    877873        domain.set_quantity('friction', 0)
    878874
    879         Br = Reflective_boundary(domain)
     875        Br = anuga.Reflective_boundary(domain)
    880876        domain.set_boundary({'exterior': Br})
    881877
     
    883879        # restricted to a polygon enclosing triangle #1 (bce)
    884880        domain.forcing_terms = []
    885         R = Rainfall(domain,
     881        R = anuga.Rainfall(domain,
    886882                     rate=lambda t: 3*t + 7,
    887883                     polygon=rainfall_poly)                     
     
    943939        domain.set_quantity('friction', 0)
    944940
    945         Br = Reflective_boundary(domain)
     941        Br = anuga.Reflective_boundary(domain)
    946942        domain.set_boundary({'exterior': Br})
    947943
     
    949945        # restricted to a polygon enclosing triangle #1 (bce)
    950946        domain.forcing_terms = []
    951         R = Rainfall(domain,
     947        R = anuga.Rainfall(domain,
    952948                     rate=lambda t: 3*t + 7,
    953949                     polygon=rainfall_poly)
     
    1000996        domain.set_quantity('friction', 0)
    1001997
    1002         Br = Reflective_boundary(domain)
     998        Br = anuga.Reflective_boundary(domain)
    1003999        domain.set_boundary({'exterior': Br})
    10041000
     
    10151011
    10161012        domain.forcing_terms = []
    1017         R = Rainfall(domain,
     1013        R = anuga.Rainfall(domain,
    10181014                     rate=main_rate,
    10191015                     polygon = [[1,1], [2,1], [2,2], [1,2]],
     
    10681064        domain.set_quantity('friction', 0)
    10691065
    1070         Br = Reflective_boundary(domain)
     1066        Br = anuga.Reflective_boundary(domain)
    10711067        domain.set_boundary({'exterior': Br})
    10721068
     
    10831079
    10841080        domain.forcing_terms = []
    1085         R = Rainfall(domain,
     1081        R = anuga.Rainfall(domain,
    10861082                     rate=main_rate,
    10871083                     polygon=[[1,1], [2,1], [2,2], [1,2]],
     
    11201116        domain.set_quantity('friction', 0)
    11211117
    1122         Br = Reflective_boundary(domain)
     1118        Br = anuga.Reflective_boundary(domain)
    11231119        domain.set_boundary({'exterior': Br})
    11241120
     
    11271123        domain.forcing_terms = []
    11281124       
    1129         I = Inflow(domain, rate=2.0, center=(1,1), radius=1)
     1125        I = anuga.Inflow(domain, rate=2.0, center=(1,1), radius=1)
    11301126        domain.forcing_terms.append(I)       
    11311127        domain.compute_forcing_terms()
     
    11541150        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    11551151
    1156         domain = Domain(points, vertices)
     1152        domain = anuga.Domain(points, vertices)
    11571153
    11581154        # Flat surface with 1m of water
     
    11611157        domain.set_quantity('friction', 0)
    11621158
    1163         Br = Reflective_boundary(domain)
     1159        Br = anuga.Reflective_boundary(domain)
    11641160        domain.set_boundary({'exterior': Br})
    11651161
     
    12081204        domain.set_quantity('friction', 0)
    12091205
    1210         Br = Reflective_boundary(domain)
     1206        Br = anuga.Reflective_boundary(domain)
    12111207        domain.set_boundary({'exterior': Br})
    12121208
     
    13521348        verbose = False
    13531349       
    1354 
    1355         #----------------------------------------------------------------------
    1356         # Import necessary modules
    1357         #----------------------------------------------------------------------
    1358 
    1359         from anuga.abstract_2d_finite_volumes.mesh_factory \
    1360                 import rectangular_cross
    1361         from anuga.shallow_water import Domain
    1362         from anuga.shallow_water.shallow_water_domain import Reflective_boundary
    1363         from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
    1364         from anuga.shallow_water.forcing import Inflow
    1365         from anuga.shallow_water.data_manager \
    1366                 import get_flow_through_cross_section
    1367         from anuga.abstract_2d_finite_volumes.util \
    1368                 import sww2csv_gauges, csv2timeseries_graphs
    1369 
    13701350        #----------------------------------------------------------------------
    13711351        # Setup computational domain
     
    14141394
    14151395                # Fixed Flowrate onto Area
    1416                 fixed_inflow = Inflow(domain,
     1396                fixed_inflow = anuga.Inflow(domain,
    14171397                                      center=(10.0, 10.0),
    14181398                                      radius=5.00,
     
    14231403                    domain.forcing_terms.append(fixed_inflow)
    14241404               
    1425                 ref_flow = fixed_inflow.rate*number_of_inflows
     1405                ref_flow = fixed_inflow.rate*number_of_inflowsg
    14261406
    14271407                # Compute normal depth on plane using Mannings equation
  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_reporting.py

    r7573 r7866  
    66import tempfile
    77
     8import anuga
     9
    810from anuga.config import g, epsilon
    911from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    1012from anuga.utilities.numerical_tools import mean
    11 from anuga.utilities.polygon import is_inside_polygon
     13from anuga.geometry.polygon import is_inside_polygon
    1214from anuga.coordinate_transforms.geo_reference import Geo_reference
    1315from anuga.abstract_2d_finite_volumes.quantity import Quantity
     
    7779        # Setup boundary conditions
    7880        #--------------------------------------------------------------
    79         Br = Reflective_boundary(domain)              # Reflective wall
     81        Br = anuga.Reflective_boundary(domain)       # Reflective wall
    8082        # Constant inflow
    81         Bd = Dirichlet_boundary([final_runup_height, 0, 0])
     83        Bd = anuga.Dirichlet_boundary([final_runup_height, 0, 0])
    8284
    8385        # All reflective to begin with (still water)
Note: See TracChangeset for help on using the changeset viewer.