Changeset 4776 for anuga_core/source/anuga/shallow_water/data_manager.py
- Timestamp:
- Nov 1, 2007, 11:12:32 AM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/data_manager.py
r4775 r4776 1290 1290 """ 1291 1291 1292 # FIXME: Can this be written feasibly using write_pts?1292 # FIXME: Can this be written feasibly using write_pts? 1293 1293 1294 1294 import os … … 1298 1298 root = basename_in 1299 1299 1300 # Get NetCDF1301 infile = NetCDFFile(root + '.dem', 'r') # Open existing netcdf file for read1300 # Get NetCDF 1301 infile = NetCDFFile(root + '.dem', 'r') # Open existing netcdf file for read 1302 1302 1303 1303 if verbose: print 'Reading DEM from %s' %(root + '.dem') … … 1305 1305 ncols = infile.ncols[0] 1306 1306 nrows = infile.nrows[0] 1307 xllcorner = infile.xllcorner[0] # Easting of lower left corner1308 yllcorner = infile.yllcorner[0] # Northing of lower left corner1307 xllcorner = infile.xllcorner[0] # Easting of lower left corner 1308 yllcorner = infile.yllcorner[0] # Northing of lower left corner 1309 1309 cellsize = infile.cellsize[0] 1310 1310 NODATA_value = infile.NODATA_value[0] … … 1315 1315 false_northing = infile.false_northing[0] 1316 1316 1317 # Text strings1317 # Text strings 1318 1318 projection = infile.projection 1319 1319 datum = infile.datum … … 1321 1321 1322 1322 1323 # Get output file1323 # Get output file 1324 1324 if basename_out == None: 1325 1325 ptsname = root + '.pts' … … 1331 1331 outfile = NetCDFFile(ptsname, 'w') 1332 1332 1333 # Create new file1333 # Create new file 1334 1334 outfile.institution = 'Geoscience Australia' 1335 1335 outfile.description = 'NetCDF pts format for compact and portable storage ' +\ 1336 1336 'of spatial point data' 1337 # assign default values1337 # Assign default values 1338 1338 if easting_min is None: easting_min = xllcorner 1339 1339 if easting_max is None: easting_max = xllcorner + ncols*cellsize … … 1341 1341 if northing_max is None: northing_max = yllcorner + nrows*cellsize 1342 1342 1343 # compute offsets to update georeferencing1343 # Compute offsets to update georeferencing 1344 1344 easting_offset = xllcorner - easting_min 1345 1345 northing_offset = yllcorner - northing_min 1346 1346 1347 # Georeferencing1347 # Georeferencing 1348 1348 outfile.zone = zone 1349 outfile.xllcorner = easting_min # Easting of lower left corner1350 outfile.yllcorner = northing_min # Northing of lower left corner1349 outfile.xllcorner = easting_min # Easting of lower left corner 1350 outfile.yllcorner = northing_min # Northing of lower left corner 1351 1351 outfile.false_easting = false_easting 1352 1352 outfile.false_northing = false_northing … … 1357 1357 1358 1358 1359 # Grid info (FIXME: probably not going to be used, but heck)1359 # Grid info (FIXME: probably not going to be used, but heck) 1360 1360 outfile.ncols = ncols 1361 1361 outfile.nrows = nrows … … 1364 1364 totalnopoints = nrows*ncols 1365 1365 1366 # calculating number of NODATA_values for each row in clipped region1367 # FIXME: use array operations to do faster1366 # Calculating number of NODATA_values for each row in clipped region 1367 # FIXME: use array operations to do faster 1368 1368 nn = 0 1369 1369 k = 0 … … 1390 1390 index2 = thisj 1391 1391 1392 # dimension definitions1392 # Dimension definitions 1393 1393 nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize)) 1394 1394 ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize)) … … 1407 1407 outfile.createDimension('number_of_dimensions', 2) #This is 2d data 1408 1408 1409 # variable definitions1409 # Variable definitions 1410 1410 outfile.createVariable('points', Float, ('number_of_points', 1411 1411 'number_of_dimensions')) … … 1417 1417 1418 1418 lenv = index2-index1+1 1419 # Store data1419 # Store data 1420 1420 global_index = 0 1421 # for i in range(nrows):1421 # for i in range(nrows): 1422 1422 for i in range(i1_0,thisi+1,1): 1423 1423 if verbose and i%((nrows+10)/10)==0: … … 1429 1429 no_NODATA = sum(v == NODATA_value) 1430 1430 if no_NODATA > 0: 1431 newcols = lenv - no_NODATA # ncols_in_bounding_box - no_NODATA1431 newcols = lenv - no_NODATA # ncols_in_bounding_box - no_NODATA 1432 1432 else: 1433 newcols = lenv # ncols_in_bounding_box1433 newcols = lenv # ncols_in_bounding_box 1434 1434 1435 1435 telev = zeros(newcols, Float)
Note: See TracChangeset
for help on using the changeset viewer.