Changeset 1866


Ignore:
Timestamp:
Oct 4, 2005, 8:48:47 PM (18 years ago)
Author:
ole
Message:

Obsoleted old versions of sww2asc and sww2ers

File:
1 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/data_manager.py

    r1865 r1866  
    1414.pts: NetCDF format for storing arbitrary points and associated attributes
    1515
    16 .asc: ASCII foramt of regular DEMs as output from ArcView
     16.asc: ASCII format of regular DEMs as output from ArcView
    1717.prj: Associated ArcView file giving more meta data for asc format
     18.ers: ERMapper header format of regular DEMs for ArcView
    1819
    1920.dem: NetCDF representation of regular DEM data
     
    12021203    infile.close()
    12031204    outfile.close()
    1204 
    1205 
    1206 def sww2asc(basename_in, basename_out = None,
    1207             quantity = None,
    1208             timestep = None,
    1209             reduction = None,
    1210             cellsize = 10,
    1211             verbose = False,
    1212             origin = None):
    1213     """Read SWW file and convert to Digitial Elevation model format (.asc)
    1214 
    1215     Example:
    1216 
    1217     ncols         3121
    1218     nrows         1800
    1219     xllcorner     722000
    1220     yllcorner     5893000
    1221     cellsize      25
    1222     NODATA_value  -9999
    1223     138.3698 137.4194 136.5062 135.5558 ..........
    1224 
    1225     Also write accompanying file with same basename_in but extension .prj
    1226     used to fix the UTM zone, datum, false northings and eastings.
    1227 
    1228     The prj format is assumed to be as
    1229 
    1230     Projection    UTM
    1231     Zone          56
    1232     Datum         WGS84
    1233     Zunits        NO
    1234     Units         METERS
    1235     Spheroid      WGS84
    1236     Xshift        0.0000000000
    1237     Yshift        10000000.0000000000
    1238     Parameters
    1239 
    1240 
    1241     if quantity is given, out values from quantity otherwise default to
    1242     elevation
    1243 
    1244     if timestep (an index) is given, output quantity at that timestep
    1245 
    1246     if reduction is given use that to reduce quantity over all timesteps.
    1247 
    1248     """
    1249     from Numeric import array, Float, concatenate, NewAxis, zeros,\
    1250          sometrue
    1251 
    1252 
    1253     #FIXME: Should be variable
    1254     datum = 'WGS84'
    1255     false_easting = 500000
    1256     false_northing = 10000000
    1257 
    1258     if quantity is None:
    1259         quantity = 'elevation'
    1260 
    1261     if reduction is None:
    1262         reduction = max
    1263 
    1264     if basename_out is None:
    1265         basename_out = basename_in + '_%s' %quantity
    1266 
    1267     swwfile = basename_in + '.sww'
    1268     ascfile = basename_out + '.asc'
    1269     prjfile = basename_out + '.prj'
    1270 
    1271 
    1272     if verbose: print 'Reading from %s' %swwfile
    1273     #Read sww file
    1274     from Scientific.IO.NetCDF import NetCDFFile
    1275     fid = NetCDFFile(swwfile)
    1276 
    1277     #Get extent and reference
    1278     x = fid.variables['x'][:]
    1279     y = fid.variables['y'][:]
    1280     volumes = fid.variables['volumes'][:]
    1281 
    1282     ymin = min(y); ymax = max(y)
    1283     xmin = min(x); xmax = max(x)
    1284 
    1285     number_of_timesteps = fid.dimensions['number_of_timesteps']
    1286     number_of_points = fid.dimensions['number_of_points']
    1287     if origin is None:
    1288 
    1289         #Get geo_reference
    1290         #sww files don't have to have a geo_ref
    1291         try:
    1292             geo_reference = Geo_reference(NetCDFObject=fid)
    1293         except AttributeError, e:
    1294             geo_reference = Geo_reference() #Default georef object
    1295 
    1296         xllcorner = geo_reference.get_xllcorner()
    1297         yllcorner = geo_reference.get_yllcorner()
    1298         zone = geo_reference.get_zone()
    1299     else:
    1300         zone = origin[0]
    1301         xllcorner = origin[1]
    1302         yllcorner = origin[2]
    1303 
    1304 
    1305     #Get quantity and reduce if applicable
    1306     if verbose: print 'Reading quantity %s' %quantity
    1307 
    1308     if quantity.lower() == 'depth':
    1309         q = fid.variables['stage'][:] - fid.variables['elevation'][:]
    1310     else:
    1311         q = fid.variables[quantity][:]
    1312 
    1313 
    1314     if len(q.shape) == 2:
    1315         if verbose: print 'Reducing quantity %s' %quantity
    1316         q_reduced = zeros( number_of_points, Float )
    1317 
    1318         for k in range(number_of_points):
    1319             q_reduced[k] = reduction( q[:,k] )
    1320 
    1321         q = q_reduced
    1322 
    1323     #Now q has dimension: number_of_points
    1324 
    1325     #Create grid and update xll/yll corner and x,y
    1326     if verbose: print 'Creating grid'
    1327     ncols = int((xmax-xmin)/cellsize)+1
    1328     nrows = int((ymax-ymin)/cellsize)+1
    1329 
    1330     newxllcorner = xmin+xllcorner
    1331     newyllcorner = ymin+yllcorner
    1332 
    1333     x = x+xllcorner-newxllcorner
    1334     y = y+yllcorner-newyllcorner
    1335 
    1336     vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
    1337     assert len(vertex_points.shape) == 2
    1338 
    1339 
    1340     from Numeric import zeros, Float
    1341     grid_points = zeros ( (ncols*nrows, 2), Float )
    1342 
    1343 
    1344     for i in xrange(nrows):
    1345         yg = i*cellsize
    1346         for j in xrange(ncols):
    1347             xg = j*cellsize
    1348             k = i*ncols + j
    1349 
    1350             grid_points[k,0] = xg
    1351             grid_points[k,1] = yg
    1352 
    1353     #Interpolate
    1354     from least_squares import Interpolation
    1355     from util import inside_polygon
    1356 
    1357     #FIXME: This should be done with precrop = True, otherwise it'll
    1358     #take forever. With expand_search  set to False, some grid points might
    1359     #miss out....
    1360 
    1361     interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
    1362                            precrop = False, expand_search = False,
    1363                            verbose = verbose)
    1364 
    1365     #Interpolate using quantity values
    1366     if verbose: print 'Interpolating'
    1367     grid_values = interp.interpolate(q).flat
    1368 
    1369     #Write
    1370     #Write prj file
    1371     if verbose: print 'Writing %s' %prjfile
    1372     prjid = open(prjfile, 'w')
    1373     prjid.write('Projection    %s\n' %'UTM')
    1374     prjid.write('Zone          %d\n' %zone)
    1375     prjid.write('Datum         %s\n' %datum)
    1376     prjid.write('Zunits        NO\n')
    1377     prjid.write('Units         METERS\n')
    1378     prjid.write('Spheroid      %s\n' %datum)
    1379     prjid.write('Xshift        %d\n' %false_easting)
    1380     prjid.write('Yshift        %d\n' %false_northing)
    1381     prjid.write('Parameters\n')
    1382     prjid.close()
    1383 
    1384 
    1385    
    1386     if verbose: print 'Writing %s' %ascfile
    1387     NODATA_value = -9999
    1388  
    1389     ascid = open(ascfile, 'w')
    1390 
    1391     ascid.write('ncols         %d\n' %ncols)
    1392     ascid.write('nrows         %d\n' %nrows)
    1393     ascid.write('xllcorner     %d\n' %newxllcorner)
    1394     ascid.write('yllcorner     %d\n' %newyllcorner)
    1395     ascid.write('cellsize      %f\n' %cellsize)
    1396     ascid.write('NODATA_value  %d\n' %NODATA_value)
    1397 
    1398 
    1399     #Get bounding polygon from mesh
    1400     P = interp.mesh.get_boundary_polygon()
    1401     inside_indices = inside_polygon(grid_points, P)
    1402 
    1403     for i in range(nrows):
    1404         if verbose and i%((nrows+10)/10)==0:
    1405             print 'Doing row %d of %d' %(i, nrows)
    1406 
    1407         for j in range(ncols):
    1408             index = (nrows-i-1)*ncols+j
    1409 
    1410             if sometrue(inside_indices == index):
    1411                 ascid.write('%f ' %grid_values[index])
    1412             else:
    1413                 ascid.write('%d ' %NODATA_value)
    1414 
    1415         ascid.write('\n')
    1416 
    1417     #Close
    1418     ascid.close()
    1419     fid.close()
    1420 
    1421 def sww2ers(basename_in, basename_out = None,
    1422             quantity = None,
    1423             timestep = None,
    1424             reduction = None,
    1425             cellsize = 10,
    1426             verbose = False,
    1427             origin = None,
    1428             datum = 'WGS84'):
    1429    
    1430     """Read SWW file and convert to Digitial Elevation model format (.asc)
    1431 
    1432     Example:
    1433 
    1434     ncols         3121
    1435     nrows         1800
    1436     xllcorner     722000
    1437     yllcorner     5893000
    1438     cellsize      25
    1439     NODATA_value  -9999
    1440     138.3698 137.4194 136.5062 135.5558 ..........
    1441 
    1442     Also write accompanying file with same basename_in but extension .prj
    1443     used to fix the UTM zone, datum, false northings and eastings.
    1444 
    1445     The prj format is assumed to be as
    1446 
    1447     Projection    UTM
    1448     Zone          56
    1449     Datum         WGS84
    1450     Zunits        NO
    1451     Units         METERS
    1452     Spheroid      WGS84
    1453     Xshift        0.0000000000
    1454     Yshift        10000000.0000000000
    1455     Parameters
    1456 
    1457 
    1458     if quantity is given, out values from quantity otherwise default to
    1459     elevation
    1460 
    1461     if timestep (an index) is given, output quantity at that timestep
    1462 
    1463     if reduction is given use that to reduce quantity over all timesteps.
    1464 
    1465     """
    1466     from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
    1467     import ermapper_grids
    1468 
    1469     header = {}
    1470     false_easting = 500000
    1471     false_northing = 10000000
    1472     NODATA_value = 0
    1473 
    1474     if quantity is None:
    1475         quantity = 'elevation'
    1476 
    1477     if reduction is None:
    1478         reduction = max
    1479 
    1480     if basename_out is None:
    1481         basename_out = basename_in + '_%s' %quantity
    1482  
    1483     swwfile = basename_in + '.sww'
    1484     # Note the use of a .ers extension is optional (write_ermapper_grid will
    1485     # deal with either option
    1486     ersfile = basename_out
    1487 ##    prjfile = basename_out + '.prj'
    1488 
    1489 
    1490     if verbose: print 'Reading from %s' %swwfile
    1491     #Read sww file
    1492     from Scientific.IO.NetCDF import NetCDFFile
    1493     fid = NetCDFFile(swwfile)
    1494 
    1495     #Get extent and reference
    1496     x = fid.variables['x'][:]
    1497     y = fid.variables['y'][:]
    1498     volumes = fid.variables['volumes'][:]
    1499 
    1500     ymin = min(y); ymax = max(y)
    1501     xmin = min(x); xmax = max(x)
    1502 
    1503     number_of_timesteps = fid.dimensions['number_of_timesteps']
    1504     number_of_points = fid.dimensions['number_of_points']
    1505     if origin is None:
    1506 
    1507         #Get geo_reference
    1508         #sww files don't have to have a geo_ref
    1509         try:
    1510             geo_reference = Geo_reference(NetCDFObject=fid)
    1511         except AttributeError, e:
    1512             geo_reference = Geo_reference() #Default georef object
    1513 
    1514         xllcorner = geo_reference.get_xllcorner()
    1515         yllcorner = geo_reference.get_yllcorner()
    1516         zone = geo_reference.get_zone()
    1517     else:
    1518         zone = origin[0]
    1519         xllcorner = origin[1]
    1520         yllcorner = origin[2]
    1521 
    1522 
    1523     #Get quantity and reduce if applicable
    1524     if verbose: print 'Reading quantity %s' %quantity
    1525 
    1526     if quantity.lower() == 'depth':
    1527         q = fid.variables['stage'][:] - fid.variables['elevation'][:]
    1528     else:
    1529         q = fid.variables[quantity][:]
    1530 
    1531 
    1532     if len(q.shape) == 2:
    1533         if verbose: print 'Reducing quantity %s' %quantity
    1534         q_reduced = zeros( number_of_points, Float )
    1535 
    1536         for k in range(number_of_points):
    1537             q_reduced[k] = reduction( q[:,k] )
    1538 
    1539         q = q_reduced
    1540 
    1541     #Now q has dimension: number_of_points
    1542 
    1543     #Create grid and update xll/yll corner and x,y
    1544     if verbose: print 'Creating grid'
    1545     ncols = int((xmax-xmin)/cellsize)+1
    1546     nrows = int((ymax-ymin)/cellsize)+1
    1547 
    1548     newxllcorner = xmin+xllcorner
    1549     newyllcorner = ymin+yllcorner
    1550 
    1551     x = x+xllcorner-newxllcorner
    1552     y = y+yllcorner-newyllcorner
    1553 
    1554     vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
    1555     assert len(vertex_points.shape) == 2
    1556 
    1557 
    1558     from Numeric import zeros, Float
    1559     grid_points = zeros ( (ncols*nrows, 2), Float )
    1560 
    1561 
    1562     for i in xrange(nrows):
    1563         yg = (nrows-i)*cellsize # this will flip the order of the y values
    1564         for j in xrange(ncols):
    1565             xg = j*cellsize
    1566             k = i*ncols + j
    1567 
    1568             grid_points[k,0] = xg
    1569             grid_points[k,1] = yg
    1570 
    1571     #Interpolate
    1572     from least_squares import Interpolation
    1573     from util import inside_polygon
    1574 
    1575     #FIXME: This should be done with precrop = True (?), otherwise it'll
    1576     #take forever. With expand_search  set to False, some grid points might
    1577     #miss out....
    1578 
    1579     interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
    1580                            precrop = False, expand_search = False,
    1581                            verbose = verbose)
    1582 
    1583     #Interpolate using quantity values
    1584     if verbose: print 'Interpolating'
    1585     grid_values = interp.interpolate(q).flat
    1586     grid_values = reshape(grid_values,(nrows, ncols))
    1587 
    1588 
    1589     # setup header information
    1590     header['datum'] = '"' + datum + '"'
    1591     # FIXME The use of hardwired UTM and zone number needs to be made optional
    1592     # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
    1593     header['projection'] = '"UTM-' + str(zone) + '"'
    1594     header['coordinatetype'] = 'EN'
    1595     if header['coordinatetype'] == 'LL':
    1596         header['longitude'] = str(newxllcorner)
    1597         header['latitude'] = str(newyllcorner)
    1598     elif header['coordinatetype'] == 'EN':
    1599         header['eastings'] = str(newxllcorner)
    1600         header['northings'] = str(newyllcorner)
    1601     header['nullcellvalue'] = str(NODATA_value)
    1602     header['xdimension'] = str(cellsize)
    1603     header['ydimension'] = str(cellsize)
    1604     header['value'] = '"' + quantity + '"'
    1605 
    1606 
    1607     #Write
    1608     if verbose: print 'Writing %s' %ersfile
    1609     ermapper_grids.write_ermapper_grid(ersfile,grid_values, header)
    1610 
    1611     fid.close()   
    1612 
    1613 ##    ascid = open(ascfile, 'w')
    1614 ##
    1615 ##    ascid.write('ncols         %d\n' %ncols)
    1616 ##    ascid.write('nrows         %d\n' %nrows)
    1617 ##    ascid.write('xllcorner     %d\n' %newxllcorner)
    1618 ##    ascid.write('yllcorner     %d\n' %newyllcorner)
    1619 ##    ascid.write('cellsize      %f\n' %cellsize)
    1620 ##    ascid.write('NODATA_value  %d\n' %NODATA_value)
    1621 ##
    1622 ##
    1623 ##    #Get bounding polygon from mesh
    1624 ##    P = interp.mesh.get_boundary_polygon()
    1625 ##    inside_indices = inside_polygon(grid_points, P)
    1626 ##
    1627 ##    for i in range(nrows):
    1628 ##        if verbose and i%((nrows+10)/10)==0:
    1629 ##            print 'Doing row %d of %d' %(i, nrows)
    1630 ##
    1631 ##        for j in range(ncols):
    1632 ##            index = (nrows-i-1)*ncols+j
    1633 ##
    1634 ##            if sometrue(inside_indices == index):
    1635 ##                ascid.write('%f ' %grid_values[index])
    1636 ##            else:
    1637 ##                ascid.write('%d ' %NODATA_value)
    1638 ##
    1639 ##        ascid.write('\n')
    1640 ##
    1641 ##    #Close
    1642 ##    ascid.close()
    1643 ##    fid.close()
    16441205
    16451206
     
    19031464        ascid.close()
    19041465        fid.close()
     1466
     1467#Backwards compatibility       
     1468def sww2asc(basename_in, basename_out = None,
     1469            quantity = None,
     1470            timestep = None,
     1471            reduction = None,
     1472            cellsize = 10,
     1473            verbose = False,
     1474            origin = None):
     1475    print 'sww2asc will soon be obsoleted - please use sww2dem'
     1476    sww2dem(basename_in,
     1477            basename_out,
     1478            quantity,
     1479            timestep,
     1480            reduction,
     1481            cellsize,
     1482            verbose,
     1483            origin,
     1484            datum = 'WGS84',
     1485            format = 'asc')
     1486           
     1487def sww2ers(basename_in, basename_out = None,
     1488            quantity = None,
     1489            timestep = None,
     1490            reduction = None,
     1491            cellsize = 10,
     1492            verbose = False,
     1493            origin = None,
     1494            datum = 'WGS84'):
     1495    print 'sww2ers will soon be obsoleted - please use sww2dem'
     1496    sww2dem(basename_in,
     1497            basename_out,
     1498            quantity,
     1499            timestep,
     1500            reduction,
     1501            cellsize,
     1502            verbose,
     1503            origin,
     1504            datum,
     1505            format = 'ers')         
     1506################################# END COMPATIBILITY ##############         
    19051507
    19061508
     
    31992801    outfile.close()
    32002802
     2803
     2804   
     2805def sww2asc_obsolete(basename_in, basename_out = None,
     2806            quantity = None,
     2807            timestep = None,
     2808            reduction = None,
     2809            cellsize = 10,
     2810            verbose = False,
     2811            origin = None):
     2812    """Read SWW file and convert to Digitial Elevation model format (.asc)
     2813
     2814    Example:
     2815
     2816    ncols         3121
     2817    nrows         1800
     2818    xllcorner     722000
     2819    yllcorner     5893000
     2820    cellsize      25
     2821    NODATA_value  -9999
     2822    138.3698 137.4194 136.5062 135.5558 ..........
     2823
     2824    Also write accompanying file with same basename_in but extension .prj
     2825    used to fix the UTM zone, datum, false northings and eastings.
     2826
     2827    The prj format is assumed to be as
     2828
     2829    Projection    UTM
     2830    Zone          56
     2831    Datum         WGS84
     2832    Zunits        NO
     2833    Units         METERS
     2834    Spheroid      WGS84
     2835    Xshift        0.0000000000
     2836    Yshift        10000000.0000000000
     2837    Parameters
     2838
     2839
     2840    if quantity is given, out values from quantity otherwise default to
     2841    elevation
     2842
     2843    if timestep (an index) is given, output quantity at that timestep
     2844
     2845    if reduction is given use that to reduce quantity over all timesteps.
     2846
     2847    """
     2848    from Numeric import array, Float, concatenate, NewAxis, zeros,\
     2849         sometrue
     2850
     2851
     2852    #FIXME: Should be variable
     2853    datum = 'WGS84'
     2854    false_easting = 500000
     2855    false_northing = 10000000
     2856
     2857    if quantity is None:
     2858        quantity = 'elevation'
     2859
     2860    if reduction is None:
     2861        reduction = max
     2862
     2863    if basename_out is None:
     2864        basename_out = basename_in + '_%s' %quantity
     2865
     2866    swwfile = basename_in + '.sww'
     2867    ascfile = basename_out + '.asc'
     2868    prjfile = basename_out + '.prj'
     2869
     2870
     2871    if verbose: print 'Reading from %s' %swwfile
     2872    #Read sww file
     2873    from Scientific.IO.NetCDF import NetCDFFile
     2874    fid = NetCDFFile(swwfile)
     2875
     2876    #Get extent and reference
     2877    x = fid.variables['x'][:]
     2878    y = fid.variables['y'][:]
     2879    volumes = fid.variables['volumes'][:]
     2880
     2881    ymin = min(y); ymax = max(y)
     2882    xmin = min(x); xmax = max(x)
     2883
     2884    number_of_timesteps = fid.dimensions['number_of_timesteps']
     2885    number_of_points = fid.dimensions['number_of_points']
     2886    if origin is None:
     2887
     2888        #Get geo_reference
     2889        #sww files don't have to have a geo_ref
     2890        try:
     2891            geo_reference = Geo_reference(NetCDFObject=fid)
     2892        except AttributeError, e:
     2893            geo_reference = Geo_reference() #Default georef object
     2894
     2895        xllcorner = geo_reference.get_xllcorner()
     2896        yllcorner = geo_reference.get_yllcorner()
     2897        zone = geo_reference.get_zone()
     2898    else:
     2899        zone = origin[0]
     2900        xllcorner = origin[1]
     2901        yllcorner = origin[2]
     2902
     2903
     2904    #Get quantity and reduce if applicable
     2905    if verbose: print 'Reading quantity %s' %quantity
     2906
     2907    if quantity.lower() == 'depth':
     2908        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
     2909    else:
     2910        q = fid.variables[quantity][:]
     2911
     2912
     2913    if len(q.shape) == 2:
     2914        if verbose: print 'Reducing quantity %s' %quantity
     2915        q_reduced = zeros( number_of_points, Float )
     2916
     2917        for k in range(number_of_points):
     2918            q_reduced[k] = reduction( q[:,k] )
     2919
     2920        q = q_reduced
     2921
     2922    #Now q has dimension: number_of_points
     2923
     2924    #Create grid and update xll/yll corner and x,y
     2925    if verbose: print 'Creating grid'
     2926    ncols = int((xmax-xmin)/cellsize)+1
     2927    nrows = int((ymax-ymin)/cellsize)+1
     2928
     2929    newxllcorner = xmin+xllcorner
     2930    newyllcorner = ymin+yllcorner
     2931
     2932    x = x+xllcorner-newxllcorner
     2933    y = y+yllcorner-newyllcorner
     2934
     2935    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
     2936    assert len(vertex_points.shape) == 2
     2937
     2938
     2939    from Numeric import zeros, Float
     2940    grid_points = zeros ( (ncols*nrows, 2), Float )
     2941
     2942
     2943    for i in xrange(nrows):
     2944        yg = i*cellsize
     2945        for j in xrange(ncols):
     2946            xg = j*cellsize
     2947            k = i*ncols + j
     2948
     2949            grid_points[k,0] = xg
     2950            grid_points[k,1] = yg
     2951
     2952    #Interpolate
     2953    from least_squares import Interpolation
     2954    from util import inside_polygon
     2955
     2956    #FIXME: This should be done with precrop = True, otherwise it'll
     2957    #take forever. With expand_search  set to False, some grid points might
     2958    #miss out....
     2959
     2960    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
     2961                           precrop = False, expand_search = False,
     2962                           verbose = verbose)
     2963
     2964    #Interpolate using quantity values
     2965    if verbose: print 'Interpolating'
     2966    grid_values = interp.interpolate(q).flat
     2967
     2968    #Write
     2969    #Write prj file
     2970    if verbose: print 'Writing %s' %prjfile
     2971    prjid = open(prjfile, 'w')
     2972    prjid.write('Projection    %s\n' %'UTM')
     2973    prjid.write('Zone          %d\n' %zone)
     2974    prjid.write('Datum         %s\n' %datum)
     2975    prjid.write('Zunits        NO\n')
     2976    prjid.write('Units         METERS\n')
     2977    prjid.write('Spheroid      %s\n' %datum)
     2978    prjid.write('Xshift        %d\n' %false_easting)
     2979    prjid.write('Yshift        %d\n' %false_northing)
     2980    prjid.write('Parameters\n')
     2981    prjid.close()
     2982
     2983
     2984   
     2985    if verbose: print 'Writing %s' %ascfile
     2986    NODATA_value = -9999
     2987 
     2988    ascid = open(ascfile, 'w')
     2989
     2990    ascid.write('ncols         %d\n' %ncols)
     2991    ascid.write('nrows         %d\n' %nrows)
     2992    ascid.write('xllcorner     %d\n' %newxllcorner)
     2993    ascid.write('yllcorner     %d\n' %newyllcorner)
     2994    ascid.write('cellsize      %f\n' %cellsize)
     2995    ascid.write('NODATA_value  %d\n' %NODATA_value)
     2996
     2997
     2998    #Get bounding polygon from mesh
     2999    P = interp.mesh.get_boundary_polygon()
     3000    inside_indices = inside_polygon(grid_points, P)
     3001
     3002    for i in range(nrows):
     3003        if verbose and i%((nrows+10)/10)==0:
     3004            print 'Doing row %d of %d' %(i, nrows)
     3005
     3006        for j in range(ncols):
     3007            index = (nrows-i-1)*ncols+j
     3008
     3009            if sometrue(inside_indices == index):
     3010                ascid.write('%f ' %grid_values[index])
     3011            else:
     3012                ascid.write('%d ' %NODATA_value)
     3013
     3014        ascid.write('\n')
     3015
     3016    #Close
     3017    ascid.close()
     3018    fid.close()
     3019
     3020def sww2ers_obsolete(basename_in, basename_out = None,
     3021            quantity = None,
     3022            timestep = None,
     3023            reduction = None,
     3024            cellsize = 10,
     3025            verbose = False,
     3026            origin = None,
     3027            datum = 'WGS84'):
     3028   
     3029    """Read SWW file and convert to Digitial Elevation model format (.asc)
     3030
     3031    Example:
     3032
     3033    ncols         3121
     3034    nrows         1800
     3035    xllcorner     722000
     3036    yllcorner     5893000
     3037    cellsize      25
     3038    NODATA_value  -9999
     3039    138.3698 137.4194 136.5062 135.5558 ..........
     3040
     3041    Also write accompanying file with same basename_in but extension .prj
     3042    used to fix the UTM zone, datum, false northings and eastings.
     3043
     3044    The prj format is assumed to be as
     3045
     3046    Projection    UTM
     3047    Zone          56
     3048    Datum         WGS84
     3049    Zunits        NO
     3050    Units         METERS
     3051    Spheroid      WGS84
     3052    Xshift        0.0000000000
     3053    Yshift        10000000.0000000000
     3054    Parameters
     3055
     3056
     3057    if quantity is given, out values from quantity otherwise default to
     3058    elevation
     3059
     3060    if timestep (an index) is given, output quantity at that timestep
     3061
     3062    if reduction is given use that to reduce quantity over all timesteps.
     3063
     3064    """
     3065    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
     3066    import ermapper_grids
     3067
     3068    header = {}
     3069    false_easting = 500000
     3070    false_northing = 10000000
     3071    NODATA_value = 0
     3072
     3073    if quantity is None:
     3074        quantity = 'elevation'
     3075
     3076    if reduction is None:
     3077        reduction = max
     3078
     3079    if basename_out is None:
     3080        basename_out = basename_in + '_%s' %quantity
     3081 
     3082    swwfile = basename_in + '.sww'
     3083    # Note the use of a .ers extension is optional (write_ermapper_grid will
     3084    # deal with either option
     3085    ersfile = basename_out
     3086##    prjfile = basename_out + '.prj'
     3087
     3088
     3089    if verbose: print 'Reading from %s' %swwfile
     3090    #Read sww file
     3091    from Scientific.IO.NetCDF import NetCDFFile
     3092    fid = NetCDFFile(swwfile)
     3093
     3094    #Get extent and reference
     3095    x = fid.variables['x'][:]
     3096    y = fid.variables['y'][:]
     3097    volumes = fid.variables['volumes'][:]
     3098
     3099    ymin = min(y); ymax = max(y)
     3100    xmin = min(x); xmax = max(x)
     3101
     3102    number_of_timesteps = fid.dimensions['number_of_timesteps']
     3103    number_of_points = fid.dimensions['number_of_points']
     3104    if origin is None:
     3105
     3106        #Get geo_reference
     3107        #sww files don't have to have a geo_ref
     3108        try:
     3109            geo_reference = Geo_reference(NetCDFObject=fid)
     3110        except AttributeError, e:
     3111            geo_reference = Geo_reference() #Default georef object
     3112
     3113        xllcorner = geo_reference.get_xllcorner()
     3114        yllcorner = geo_reference.get_yllcorner()
     3115        zone = geo_reference.get_zone()
     3116    else:
     3117        zone = origin[0]
     3118        xllcorner = origin[1]
     3119        yllcorner = origin[2]
     3120
     3121
     3122    #Get quantity and reduce if applicable
     3123    if verbose: print 'Reading quantity %s' %quantity
     3124
     3125    if quantity.lower() == 'depth':
     3126        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
     3127    else:
     3128        q = fid.variables[quantity][:]
     3129
     3130
     3131    if len(q.shape) == 2:
     3132        if verbose: print 'Reducing quantity %s' %quantity
     3133        q_reduced = zeros( number_of_points, Float )
     3134
     3135        for k in range(number_of_points):
     3136            q_reduced[k] = reduction( q[:,k] )
     3137
     3138        q = q_reduced
     3139
     3140    #Now q has dimension: number_of_points
     3141
     3142    #Create grid and update xll/yll corner and x,y
     3143    if verbose: print 'Creating grid'
     3144    ncols = int((xmax-xmin)/cellsize)+1
     3145    nrows = int((ymax-ymin)/cellsize)+1
     3146
     3147    newxllcorner = xmin+xllcorner
     3148    newyllcorner = ymin+yllcorner
     3149
     3150    x = x+xllcorner-newxllcorner
     3151    y = y+yllcorner-newyllcorner
     3152
     3153    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
     3154    assert len(vertex_points.shape) == 2
     3155
     3156
     3157    from Numeric import zeros, Float
     3158    grid_points = zeros ( (ncols*nrows, 2), Float )
     3159
     3160
     3161    for i in xrange(nrows):
     3162        yg = (nrows-i)*cellsize # this will flip the order of the y values
     3163        for j in xrange(ncols):
     3164            xg = j*cellsize
     3165            k = i*ncols + j
     3166
     3167            grid_points[k,0] = xg
     3168            grid_points[k,1] = yg
     3169
     3170    #Interpolate
     3171    from least_squares import Interpolation
     3172    from util import inside_polygon
     3173
     3174    #FIXME: This should be done with precrop = True (?), otherwise it'll
     3175    #take forever. With expand_search  set to False, some grid points might
     3176    #miss out....
     3177
     3178    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
     3179                           precrop = False, expand_search = False,
     3180                           verbose = verbose)
     3181
     3182    #Interpolate using quantity values
     3183    if verbose: print 'Interpolating'
     3184    grid_values = interp.interpolate(q).flat
     3185    grid_values = reshape(grid_values,(nrows, ncols))
     3186
     3187
     3188    # setup header information
     3189    header['datum'] = '"' + datum + '"'
     3190    # FIXME The use of hardwired UTM and zone number needs to be made optional
     3191    # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
     3192    header['projection'] = '"UTM-' + str(zone) + '"'
     3193    header['coordinatetype'] = 'EN'
     3194    if header['coordinatetype'] == 'LL':
     3195        header['longitude'] = str(newxllcorner)
     3196        header['latitude'] = str(newyllcorner)
     3197    elif header['coordinatetype'] == 'EN':
     3198        header['eastings'] = str(newxllcorner)
     3199        header['northings'] = str(newyllcorner)
     3200    header['nullcellvalue'] = str(NODATA_value)
     3201    header['xdimension'] = str(cellsize)
     3202    header['ydimension'] = str(cellsize)
     3203    header['value'] = '"' + quantity + '"'
     3204
     3205
     3206    #Write
     3207    if verbose: print 'Writing %s' %ersfile
     3208    ermapper_grids.write_ermapper_grid(ersfile,grid_values, header)
     3209
     3210    fid.close()   
     3211
     3212##    ascid = open(ascfile, 'w')
     3213##
     3214##    ascid.write('ncols         %d\n' %ncols)
     3215##    ascid.write('nrows         %d\n' %nrows)
     3216##    ascid.write('xllcorner     %d\n' %newxllcorner)
     3217##    ascid.write('yllcorner     %d\n' %newyllcorner)
     3218##    ascid.write('cellsize      %f\n' %cellsize)
     3219##    ascid.write('NODATA_value  %d\n' %NODATA_value)
     3220##
     3221##
     3222##    #Get bounding polygon from mesh
     3223##    P = interp.mesh.get_boundary_polygon()
     3224##    inside_indices = inside_polygon(grid_points, P)
     3225##
     3226##    for i in range(nrows):
     3227##        if verbose and i%((nrows+10)/10)==0:
     3228##            print 'Doing row %d of %d' %(i, nrows)
     3229##
     3230##        for j in range(ncols):
     3231##            index = (nrows-i-1)*ncols+j
     3232##
     3233##            if sometrue(inside_indices == index):
     3234##                ascid.write('%f ' %grid_values[index])
     3235##            else:
     3236##                ascid.write('%d ' %NODATA_value)
     3237##
     3238##        ascid.write('\n')
     3239##
     3240##    #Close
     3241##    ascid.close()
     3242##    fid.close()
     3243
Note: See TracChangeset for help on using the changeset viewer.