Ignore:
Timestamp:
May 21, 2007, 10:15:58 AM (17 years ago)
Author:
duncan
Message:

adding export_grid

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r4453 r4462  
    11051105
    11061106
    1107 
    1108     def test_sww2dem_asc_elevation(self):
     1107    def test_sww2dem_asc_elevation_depth(self):
    11091108        """Test that sww information can be converted correctly to asc/prj
    11101109        format readable by e.g. ArcView
    … …  
    11261125        self.domain.smooth = True
    11271126        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     1127        self.domain.set_quantity('stage', 1.0)
    11281128
    11291129        self.domain.geo_reference = Geo_reference(56,308500,6189000)
    … …  
    11511151        stage = fid.variables['stage'][:]
    11521152
     1153        fid.close()
    11531154
    11541155        #Export to ascii/prj files
    … …  
    12351236            for i in range(5):
    12361237                assert allclose(float(L[i]), -i*cellsize - y)
    1237 
    1238 
    1239         fid.close()
     1238               
     1239        #Cleanup
     1240        os.remove(prjfile)
     1241        os.remove(ascfile)
     1242
     1243        #Export to ascii/prj files
     1244        sww2dem(self.domain.get_name(),
     1245                quantity = 'depth',
     1246                cellsize = cellsize,
     1247                verbose = self.verbose,
     1248                format = 'asc')
     1249       
     1250        #Check asc file
     1251        ascfile = self.domain.get_name() + '_depth.asc'
     1252        prjfile = self.domain.get_name() + '_depth.prj'
     1253        ascid = open(ascfile)
     1254        lines = ascid.readlines()
     1255        ascid.close()
     1256
     1257        L = lines[0].strip().split()
     1258        assert L[0].strip().lower() == 'ncols'
     1259        assert L[1].strip().lower() == '5'
     1260
     1261        L = lines[1].strip().split()
     1262        assert L[0].strip().lower() == 'nrows'
     1263        assert L[1].strip().lower() == '5'
     1264
     1265        L = lines[2].strip().split()
     1266        assert L[0].strip().lower() == 'xllcorner'
     1267        assert allclose(float(L[1].strip().lower()), 308500)
     1268
     1269        L = lines[3].strip().split()
     1270        assert L[0].strip().lower() == 'yllcorner'
     1271        assert allclose(float(L[1].strip().lower()), 6189000)
     1272
     1273        L = lines[4].strip().split()
     1274        assert L[0].strip().lower() == 'cellsize'
     1275        assert allclose(float(L[1].strip().lower()), cellsize)
     1276
     1277        L = lines[5].strip().split()
     1278        assert L[0].strip() == 'NODATA_value'
     1279        assert L[1].strip().lower() == '-9999'
     1280
     1281        #Check grid values
     1282        for j in range(5):
     1283            L = lines[6+j].strip().split()
     1284            y = (4-j) * cellsize
     1285            for i in range(5):
     1286                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
     1287
    12401288
    12411289        #Cleanup
    … …  
    12451293
    12461294
     1295    def test_export_grid(self):
     1296        """Test that sww information can be converted correctly to asc/prj
     1297        format readable by e.g. ArcView
     1298        """
     1299
     1300        import time, os
     1301        from Numeric import array, zeros, allclose, Float, concatenate
     1302        from Scientific.IO.NetCDF import NetCDFFile
     1303
     1304        #Setup
     1305        self.domain.set_name('datatest')
     1306
     1307        prjfile = self.domain.get_name() + '_elevation.prj'
     1308        ascfile = self.domain.get_name() + '_elevation.asc'
     1309        swwfile = self.domain.get_name() + '.sww'
     1310
     1311        self.domain.set_datadir('.')
     1312        self.domain.format = 'sww'
     1313        self.domain.smooth = True
     1314        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     1315        self.domain.set_quantity('stage', 1.0)
     1316
     1317        self.domain.geo_reference = Geo_reference(56,308500,6189000)
     1318
     1319        sww = get_dataobject(self.domain)
     1320        sww.store_connectivity()
     1321        sww.store_timestep('stage')
     1322        self.domain.evolve_to_end(finaltime = 0.01)
     1323        sww.store_timestep('stage')
     1324
     1325        cellsize = 0.25
     1326        #Check contents
     1327        #Get NetCDF
     1328
     1329        fid = NetCDFFile(sww.filename, 'r')
     1330
     1331        # Get the variables
     1332        x = fid.variables['x'][:]
     1333        y = fid.variables['y'][:]
     1334        z = fid.variables['elevation'][:]
     1335        time = fid.variables['time'][:]
     1336        stage = fid.variables['stage'][:]
     1337
     1338        fid.close()
     1339
     1340        #Export to ascii/prj files
     1341        export_grid(self.domain.get_name(),
     1342                quantities = 'elevation',
     1343                cellsize = cellsize,
     1344                verbose = self.verbose,
     1345                format = 'asc')
     1346
     1347        #Check asc file
     1348        ascid = open(ascfile)
     1349        lines = ascid.readlines()
     1350        ascid.close()
     1351
     1352        L = lines[2].strip().split()
     1353        assert L[0].strip().lower() == 'xllcorner'
     1354        assert allclose(float(L[1].strip().lower()), 308500)
     1355
     1356        L = lines[3].strip().split()
     1357        assert L[0].strip().lower() == 'yllcorner'
     1358        assert allclose(float(L[1].strip().lower()), 6189000)
     1359
     1360        #Check grid values
     1361        for j in range(5):
     1362            L = lines[6+j].strip().split()
     1363            y = (4-j) * cellsize
     1364            for i in range(5):
     1365                assert allclose(float(L[i]), -i*cellsize - y)
     1366               
     1367        #Cleanup
     1368        os.remove(prjfile)
     1369        os.remove(ascfile)
     1370        os.remove(swwfile)
     1371
     1372    def test_export_gridII(self):
     1373        """Test that sww information can be converted correctly to asc/prj
     1374        format readable by e.g. ArcView
     1375        """
     1376
     1377        import time, os
     1378        from Numeric import array, zeros, allclose, Float, concatenate
     1379        from Scientific.IO.NetCDF import NetCDFFile
     1380
     1381        #Setup
     1382        self.domain.set_name('datatest')
     1383
     1384        swwfile = self.domain.get_name() + '.sww'
     1385
     1386        self.domain.set_datadir('.')
     1387        self.domain.format = 'sww'
     1388        self.domain.smooth = True
     1389        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     1390        self.domain.set_quantity('stage', 1.0)
     1391
     1392        self.domain.geo_reference = Geo_reference(56,308500,6189000)
     1393
     1394        sww = get_dataobject(self.domain)
     1395        sww.store_connectivity()
     1396        sww.store_timestep('stage')
     1397        self.domain.evolve_to_end(finaltime = 0.01)
     1398        sww.store_timestep('stage')
     1399
     1400        cellsize = 0.25
     1401        #Check contents
     1402        #Get NetCDF
     1403
     1404        fid = NetCDFFile(sww.filename, 'r')
     1405
     1406        # Get the variables
     1407        x = fid.variables['x'][:]
     1408        y = fid.variables['y'][:]
     1409        z = fid.variables['elevation'][:]
     1410        time = fid.variables['time'][:]
     1411        stage = fid.variables['stage'][:]
     1412
     1413        fid.close()
     1414
     1415        #Export to ascii/prj files
     1416        if True:
     1417            export_grid(self.domain.get_name(),
     1418                        quantities = ['elevation', 'depth'],
     1419                        cellsize = cellsize,
     1420                        verbose = self.verbose,
     1421                        format = 'asc')
     1422
     1423        else:
     1424            export_grid(self.domain.get_name(),
     1425                quantities = ['depth'],
     1426                cellsize = cellsize,
     1427                verbose = self.verbose,
     1428                format = 'asc')
     1429
     1430
     1431            export_grid(self.domain.get_name(),
     1432                quantities = ['elevation'],
     1433                cellsize = cellsize,
     1434                verbose = self.verbose,
     1435                format = 'asc')
     1436
     1437        prjfile = self.domain.get_name() + '_elevation.prj'
     1438        ascfile = self.domain.get_name() + '_elevation.asc'
     1439       
     1440        #Check asc file
     1441        ascid = open(ascfile)
     1442        lines = ascid.readlines()
     1443        ascid.close()
     1444
     1445        L = lines[2].strip().split()
     1446        assert L[0].strip().lower() == 'xllcorner'
     1447        assert allclose(float(L[1].strip().lower()), 308500)
     1448
     1449        L = lines[3].strip().split()
     1450        assert L[0].strip().lower() == 'yllcorner'
     1451        assert allclose(float(L[1].strip().lower()), 6189000)
     1452
     1453        #print "ascfile", ascfile
     1454        #Check grid values
     1455        for j in range(5):
     1456            L = lines[6+j].strip().split()
     1457            y = (4-j) * cellsize
     1458            for i in range(5):
     1459                #print " -i*cellsize - y",  -i*cellsize - y
     1460                #print "float(L[i])", float(L[i])
     1461                assert allclose(float(L[i]), -i*cellsize - y)
     1462               
     1463        #Cleanup
     1464        os.remove(prjfile)
     1465        os.remove(ascfile)
     1466       
     1467        #Check asc file
     1468        ascfile = self.domain.get_name() + '_depth.asc'
     1469        prjfile = self.domain.get_name() + '_depth.prj'
     1470        ascid = open(ascfile)
     1471        lines = ascid.readlines()
     1472        ascid.close()
     1473
     1474        L = lines[2].strip().split()
     1475        assert L[0].strip().lower() == 'xllcorner'
     1476        assert allclose(float(L[1].strip().lower()), 308500)
     1477
     1478        L = lines[3].strip().split()
     1479        assert L[0].strip().lower() == 'yllcorner'
     1480        assert allclose(float(L[1].strip().lower()), 6189000)
     1481
     1482        #Check grid values
     1483        for j in range(5):
     1484            L = lines[6+j].strip().split()
     1485            y = (4-j) * cellsize
     1486            for i in range(5):
     1487                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
     1488
     1489        #Cleanup
     1490        os.remove(prjfile)
     1491        os.remove(ascfile)
     1492        os.remove(swwfile)
     1493
     1494
     1495    def test_export_gridIII(self):
     1496        """Test that sww information can be converted correctly to asc/prj
     1497        format readable by e.g. ArcView
     1498        """
     1499
     1500        import time, os
     1501        from Numeric import array, zeros, allclose, Float, concatenate
     1502        from Scientific.IO.NetCDF import NetCDFFile
     1503
     1504        #Setup
     1505        self.domain.set_name('datatest')
     1506
     1507        swwfile = self.domain.get_name() + '.sww'
     1508
     1509        self.domain.set_datadir('.')
     1510        self.domain.format = 'sww'
     1511        self.domain.smooth = True
     1512        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     1513        self.domain.set_quantity('stage', 1.0)
     1514
     1515        self.domain.geo_reference = Geo_reference(56,308500,6189000)
     1516
     1517        sww = get_dataobject(self.domain)
     1518        sww.store_connectivity()
     1519        sww.store_timestep('stage')
     1520        self.domain.evolve_to_end(finaltime = 0.01)
     1521        sww.store_timestep('stage')
     1522
     1523        cellsize = 0.25
     1524        #Check contents
     1525        #Get NetCDF
     1526
     1527        fid = NetCDFFile(sww.filename, 'r')
     1528
     1529        # Get the variables
     1530        x = fid.variables['x'][:]
     1531        y = fid.variables['y'][:]
     1532        z = fid.variables['elevation'][:]
     1533        time = fid.variables['time'][:]
     1534        stage = fid.variables['stage'][:]
     1535
     1536        fid.close()
     1537
     1538        #Export to ascii/prj files
     1539        extra_name_out = 'yeah'
     1540        if True:
     1541            export_grid(self.domain.get_name(),
     1542                        quantities = ['elevation', 'depth'],
     1543                        extra_name_out = extra_name_out,
     1544                        cellsize = cellsize,
     1545                        verbose = self.verbose,
     1546                        format = 'asc')
     1547
     1548        else:
     1549            export_grid(self.domain.get_name(),
     1550                quantities = ['depth'],
     1551                cellsize = cellsize,
     1552                verbose = self.verbose,
     1553                format = 'asc')
     1554
     1555
     1556            export_grid(self.domain.get_name(),
     1557                quantities = ['elevation'],
     1558                cellsize = cellsize,
     1559                verbose = self.verbose,
     1560                format = 'asc')
     1561
     1562        prjfile = self.domain.get_name() + '_elevation_yeah.prj'
     1563        ascfile = self.domain.get_name() + '_elevation_yeah.asc'
     1564       
     1565        #Check asc file
     1566        ascid = open(ascfile)
     1567        lines = ascid.readlines()
     1568        ascid.close()
     1569
     1570        L = lines[2].strip().split()
     1571        assert L[0].strip().lower() == 'xllcorner'
     1572        assert allclose(float(L[1].strip().lower()), 308500)
     1573
     1574        L = lines[3].strip().split()
     1575        assert L[0].strip().lower() == 'yllcorner'
     1576        assert allclose(float(L[1].strip().lower()), 6189000)
     1577
     1578        #print "ascfile", ascfile
     1579        #Check grid values
     1580        for j in range(5):
     1581            L = lines[6+j].strip().split()
     1582            y = (4-j) * cellsize
     1583            for i in range(5):
     1584                #print " -i*cellsize - y",  -i*cellsize - y
     1585                #print "float(L[i])", float(L[i])
     1586                assert allclose(float(L[i]), -i*cellsize - y)
     1587               
     1588        #Cleanup
     1589        os.remove(prjfile)
     1590        os.remove(ascfile)
     1591       
     1592        #Check asc file
     1593        ascfile = self.domain.get_name() + '_depth_yeah.asc'
     1594        prjfile = self.domain.get_name() + '_depth_yeah.prj'
     1595        ascid = open(ascfile)
     1596        lines = ascid.readlines()
     1597        ascid.close()
     1598
     1599        L = lines[2].strip().split()
     1600        assert L[0].strip().lower() == 'xllcorner'
     1601        assert allclose(float(L[1].strip().lower()), 308500)
     1602
     1603        L = lines[3].strip().split()
     1604        assert L[0].strip().lower() == 'yllcorner'
     1605        assert allclose(float(L[1].strip().lower()), 6189000)
     1606
     1607        #Check grid values
     1608        for j in range(5):
     1609            L = lines[6+j].strip().split()
     1610            y = (4-j) * cellsize
     1611            for i in range(5):
     1612                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
     1613
     1614        #Cleanup
     1615        os.remove(prjfile)
     1616        os.remove(ascfile)
     1617        os.remove(swwfile)
     1618
     1619    def test_export_grid_parallel(self):
     1620        """Test that sww information can be converted correctly to asc/prj
     1621        format readable by e.g. ArcView
     1622        """
     1623
     1624        import time, os
     1625        from Numeric import array, zeros, allclose, Float, concatenate
     1626        from Scientific.IO.NetCDF import NetCDFFile
     1627
     1628        base_name = 'tegp'
     1629        #Setup
     1630        self.domain.set_name(base_name+'_P0_8')
     1631        swwfile = self.domain.get_name() + '.sww'
     1632
     1633        self.domain.set_datadir('.')
     1634        self.domain.format = 'sww'
     1635        self.domain.smooth = True
     1636        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     1637        self.domain.set_quantity('stage', 1.0)
     1638
     1639        self.domain.geo_reference = Geo_reference(56,308500,6189000)
     1640
     1641        sww = get_dataobject(self.domain)
     1642        sww.store_connectivity()
     1643        sww.store_timestep('stage')
     1644        self.domain.evolve_to_end(finaltime = 0.0001)
     1645        #Setup
     1646        self.domain.set_name(base_name+'_P1_8')
     1647        swwfile2 = self.domain.get_name() + '.sww'
     1648        sww = get_dataobject(self.domain)
     1649        sww.store_connectivity()
     1650        sww.store_timestep('stage')
     1651        self.domain.evolve_to_end(finaltime = 0.0002)
     1652        sww.store_timestep('stage')
     1653
     1654        cellsize = 0.25
     1655        #Check contents
     1656        #Get NetCDF
     1657
     1658        fid = NetCDFFile(sww.filename, 'r')
     1659
     1660        # Get the variables
     1661        x = fid.variables['x'][:]
     1662        y = fid.variables['y'][:]
     1663        z = fid.variables['elevation'][:]
     1664        time = fid.variables['time'][:]
     1665        stage = fid.variables['stage'][:]
     1666
     1667        fid.close()
     1668
     1669        #Export to ascii/prj files
     1670        extra_name_out = 'yeah'
     1671        export_grid(base_name,
     1672                    quantities = ['elevation', 'depth'],
     1673                    extra_name_out = extra_name_out,
     1674                    cellsize = cellsize,
     1675                    verbose = self.verbose,
     1676                    format = 'asc')
     1677
     1678        prjfile = base_name + '_P0_8_elevation_yeah.prj'
     1679        ascfile = base_name + '_P0_8_elevation_yeah.asc'       
     1680        #Check asc file
     1681        ascid = open(ascfile)
     1682        lines = ascid.readlines()
     1683        ascid.close()
     1684        #Check grid values
     1685        for j in range(5):
     1686            L = lines[6+j].strip().split()
     1687            y = (4-j) * cellsize
     1688            for i in range(5):
     1689                #print " -i*cellsize - y",  -i*cellsize - y
     1690                #print "float(L[i])", float(L[i])
     1691                assert allclose(float(L[i]), -i*cellsize - y)               
     1692        #Cleanup
     1693        os.remove(prjfile)
     1694        os.remove(ascfile)
     1695
     1696        prjfile = base_name + '_P1_8_elevation_yeah.prj'
     1697        ascfile = base_name + '_P1_8_elevation_yeah.asc'       
     1698        #Check asc file
     1699        ascid = open(ascfile)
     1700        lines = ascid.readlines()
     1701        ascid.close()
     1702        #Check grid values
     1703        for j in range(5):
     1704            L = lines[6+j].strip().split()
     1705            y = (4-j) * cellsize
     1706            for i in range(5):
     1707                #print " -i*cellsize - y",  -i*cellsize - y
     1708                #print "float(L[i])", float(L[i])
     1709                assert allclose(float(L[i]), -i*cellsize - y)               
     1710        #Cleanup
     1711        os.remove(prjfile)
     1712        os.remove(ascfile)
     1713        os.remove(swwfile)
     1714
     1715        #Check asc file
     1716        ascfile = base_name + '_P0_8_depth_yeah.asc'
     1717        prjfile = base_name + '_P0_8_depth_yeah.prj'
     1718        ascid = open(ascfile)
     1719        lines = ascid.readlines()
     1720        ascid.close()
     1721        #Check grid values
     1722        for j in range(5):
     1723            L = lines[6+j].strip().split()
     1724            y = (4-j) * cellsize
     1725            for i in range(5):
     1726                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
     1727        #Cleanup
     1728        os.remove(prjfile)
     1729        os.remove(ascfile)
     1730
     1731        #Check asc file
     1732        ascfile = base_name + '_P1_8_depth_yeah.asc'
     1733        prjfile = base_name + '_P1_8_depth_yeah.prj'
     1734        ascid = open(ascfile)
     1735        lines = ascid.readlines()
     1736        ascid.close()
     1737        #Check grid values
     1738        for j in range(5):
     1739            L = lines[6+j].strip().split()
     1740            y = (4-j) * cellsize
     1741            for i in range(5):
     1742                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
     1743        #Cleanup
     1744        os.remove(prjfile)
     1745        os.remove(ascfile)
     1746        os.remove(swwfile2)
    12471747
    12481748    def test_sww2dem_larger(self):
    … …  
    61876687    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sww_origin')
    61886688    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_header')
    6189     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_range')
     6689    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_grid_parallel')
     6690    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_grid')
    61906691    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
    61916692        Test_Data_Manager.verbose=True
Note: See TracChangeset for help on using the changeset viewer.