Changeset 1866
 Timestamp:
 Oct 4, 2005, 8:48:47 PM (18 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

inundation/pyvolution/data_manager.py
r1865 r1866 14 14 .pts: NetCDF format for storing arbitrary points and associated attributes 15 15 16 .asc: ASCII for amt of regular DEMs as output from ArcView16 .asc: ASCII format of regular DEMs as output from ArcView 17 17 .prj: Associated ArcView file giving more meta data for asc format 18 .ers: ERMapper header format of regular DEMs for ArcView 18 19 19 20 .dem: NetCDF representation of regular DEM data … … 1202 1203 infile.close() 1203 1204 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 31211218 nrows 18001219 xllcorner 7220001220 yllcorner 58930001221 cellsize 251222 NODATA_value 99991223 138.3698 137.4194 136.5062 135.5558 ..........1224 1225 Also write accompanying file with same basename_in but extension .prj1226 used to fix the UTM zone, datum, false northings and eastings.1227 1228 The prj format is assumed to be as1229 1230 Projection UTM1231 Zone 561232 Datum WGS841233 Zunits NO1234 Units METERS1235 Spheroid WGS841236 Xshift 0.00000000001237 Yshift 10000000.00000000001238 Parameters1239 1240 1241 if quantity is given, out values from quantity otherwise default to1242 elevation1243 1244 if timestep (an index) is given, output quantity at that timestep1245 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 sometrue1251 1252 1253 #FIXME: Should be variable1254 datum = 'WGS84'1255 false_easting = 5000001256 false_northing = 100000001257 1258 if quantity is None:1259 quantity = 'elevation'1260 1261 if reduction is None:1262 reduction = max1263 1264 if basename_out is None:1265 basename_out = basename_in + '_%s' %quantity1266 1267 swwfile = basename_in + '.sww'1268 ascfile = basename_out + '.asc'1269 prjfile = basename_out + '.prj'1270 1271 1272 if verbose: print 'Reading from %s' %swwfile1273 #Read sww file1274 from Scientific.IO.NetCDF import NetCDFFile1275 fid = NetCDFFile(swwfile)1276 1277 #Get extent and reference1278 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_reference1290 #sww files don't have to have a geo_ref1291 try:1292 geo_reference = Geo_reference(NetCDFObject=fid)1293 except AttributeError, e:1294 geo_reference = Geo_reference() #Default georef object1295 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 applicable1306 if verbose: print 'Reading quantity %s' %quantity1307 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' %quantity1316 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_reduced1322 1323 #Now q has dimension: number_of_points1324 1325 #Create grid and update xll/yll corner and x,y1326 if verbose: print 'Creating grid'1327 ncols = int((xmaxxmin)/cellsize)+11328 nrows = int((ymaxymin)/cellsize)+11329 1330 newxllcorner = xmin+xllcorner1331 newyllcorner = ymin+yllcorner1332 1333 x = x+xllcornernewxllcorner1334 y = y+yllcornernewyllcorner1335 1336 vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)1337 assert len(vertex_points.shape) == 21338 1339 1340 from Numeric import zeros, Float1341 grid_points = zeros ( (ncols*nrows, 2), Float )1342 1343 1344 for i in xrange(nrows):1345 yg = i*cellsize1346 for j in xrange(ncols):1347 xg = j*cellsize1348 k = i*ncols + j1349 1350 grid_points[k,0] = xg1351 grid_points[k,1] = yg1352 1353 #Interpolate1354 from least_squares import Interpolation1355 from util import inside_polygon1356 1357 #FIXME: This should be done with precrop = True, otherwise it'll1358 #take forever. With expand_search set to False, some grid points might1359 #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 values1366 if verbose: print 'Interpolating'1367 grid_values = interp.interpolate(q).flat1368 1369 #Write1370 #Write prj file1371 if verbose: print 'Writing %s' %prjfile1372 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' %ascfile1387 NODATA_value = 99991388 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 mesh1400 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 = (nrowsi1)*ncols+j1409 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 #Close1418 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 31211435 nrows 18001436 xllcorner 7220001437 yllcorner 58930001438 cellsize 251439 NODATA_value 99991440 138.3698 137.4194 136.5062 135.5558 ..........1441 1442 Also write accompanying file with same basename_in but extension .prj1443 used to fix the UTM zone, datum, false northings and eastings.1444 1445 The prj format is assumed to be as1446 1447 Projection UTM1448 Zone 561449 Datum WGS841450 Zunits NO1451 Units METERS1452 Spheroid WGS841453 Xshift 0.00000000001454 Yshift 10000000.00000000001455 Parameters1456 1457 1458 if quantity is given, out values from quantity otherwise default to1459 elevation1460 1461 if timestep (an index) is given, output quantity at that timestep1462 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, sometrue1467 import ermapper_grids1468 1469 header = {}1470 false_easting = 5000001471 false_northing = 100000001472 NODATA_value = 01473 1474 if quantity is None:1475 quantity = 'elevation'1476 1477 if reduction is None:1478 reduction = max1479 1480 if basename_out is None:1481 basename_out = basename_in + '_%s' %quantity1482 1483 swwfile = basename_in + '.sww'1484 # Note the use of a .ers extension is optional (write_ermapper_grid will1485 # deal with either option1486 ersfile = basename_out1487 ## prjfile = basename_out + '.prj'1488 1489 1490 if verbose: print 'Reading from %s' %swwfile1491 #Read sww file1492 from Scientific.IO.NetCDF import NetCDFFile1493 fid = NetCDFFile(swwfile)1494 1495 #Get extent and reference1496 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_reference1508 #sww files don't have to have a geo_ref1509 try:1510 geo_reference = Geo_reference(NetCDFObject=fid)1511 except AttributeError, e:1512 geo_reference = Geo_reference() #Default georef object1513 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 applicable1524 if verbose: print 'Reading quantity %s' %quantity1525 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' %quantity1534 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_reduced1540 1541 #Now q has dimension: number_of_points1542 1543 #Create grid and update xll/yll corner and x,y1544 if verbose: print 'Creating grid'1545 ncols = int((xmaxxmin)/cellsize)+11546 nrows = int((ymaxymin)/cellsize)+11547 1548 newxllcorner = xmin+xllcorner1549 newyllcorner = ymin+yllcorner1550 1551 x = x+xllcornernewxllcorner1552 y = y+yllcornernewyllcorner1553 1554 vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)1555 assert len(vertex_points.shape) == 21556 1557 1558 from Numeric import zeros, Float1559 grid_points = zeros ( (ncols*nrows, 2), Float )1560 1561 1562 for i in xrange(nrows):1563 yg = (nrowsi)*cellsize # this will flip the order of the y values1564 for j in xrange(ncols):1565 xg = j*cellsize1566 k = i*ncols + j1567 1568 grid_points[k,0] = xg1569 grid_points[k,1] = yg1570 1571 #Interpolate1572 from least_squares import Interpolation1573 from util import inside_polygon1574 1575 #FIXME: This should be done with precrop = True (?), otherwise it'll1576 #take forever. With expand_search set to False, some grid points might1577 #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 values1584 if verbose: print 'Interpolating'1585 grid_values = interp.interpolate(q).flat1586 grid_values = reshape(grid_values,(nrows, ncols))1587 1588 1589 # setup header information1590 header['datum'] = '"' + datum + '"'1591 # FIXME The use of hardwired UTM and zone number needs to be made optional1592 # 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 #Write1608 if verbose: print 'Writing %s' %ersfile1609 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 mesh1624 ## 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 = (nrowsi1)*ncols+j1633 ##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 ## #Close1642 ## ascid.close()1643 ## fid.close()1644 1205 1645 1206 … … 1903 1464 ascid.close() 1904 1465 fid.close() 1466 1467 #Backwards compatibility 1468 def 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 1487 def 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 ############## 1905 1507 1906 1508 … … 3199 2801 outfile.close() 3200 2802 2803 2804 2805 def 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((xmaxxmin)/cellsize)+1 2927 nrows = int((ymaxymin)/cellsize)+1 2928 2929 newxllcorner = xmin+xllcorner 2930 newyllcorner = ymin+yllcorner 2931 2932 x = x+xllcornernewxllcorner 2933 y = y+yllcornernewyllcorner 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 = (nrowsi1)*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 3020 def 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((xmaxxmin)/cellsize)+1 3145 nrows = int((ymaxymin)/cellsize)+1 3146 3147 newxllcorner = xmin+xllcorner 3148 newyllcorner = ymin+yllcorner 3149 3150 x = x+xllcornernewxllcorner 3151 y = y+yllcornernewyllcorner 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 = (nrowsi)*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 = (nrowsi1)*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.