Changeset 7758


Ignore:
Timestamp:
May 31, 2010, 4:45:39 PM (14 years ago)
Author:
hudson
Message:

Added extra file conversion classes to file_conversion.

Location:
trunk/anuga_core/source/anuga/file_conversion
Files:
8 added
3 edited

Legend:

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

    r7751 r7758  
    77
    88
    9 #from shallow_water_balanced_domain import Swb_domain
    10 
    11 
    12  
     9# commonly used imports
     10from anuga.file_conversion.file_conversion import sww2obj, dat2obj, \
     11                    timefile2netcdf, tsh2sww, urs2sww, urs2nc
     12                   
     13from anuga.file_conversion.dem2pts import dem2pts                   
     14from anuga.file_conversion.esri2sww import esri2sww   
     15from anuga.file_conversion.sww2dem import sww2dem     
     16from anuga.file_conversion.asc2dem import asc2dem     
     17from anuga.file_conversion.ferret2sww import ferret2sww     
  • trunk/anuga_core/source/anuga/file_conversion/file_conversion.py

    r7743 r7758  
    316316    sww.store_connectivity()
    317317    sww.store_timestep()
    318 
    319 
    320 ##
    321 # @brief Convert CSIRO ESRI file to an SWW boundary file.
    322 # @param bath_dir
    323 # @param elevation_dir
    324 # @param ucur_dir
    325 # @param vcur_dir
    326 # @param sww_file
    327 # @param minlat
    328 # @param maxlat
    329 # @param minlon
    330 # @param maxlon
    331 # @param zscale
    332 # @param mean_stage
    333 # @param fail_on_NaN
    334 # @param elevation_NaN_filler
    335 # @param bath_prefix
    336 # @param elevation_prefix
    337 # @param verbose
    338 # @note Also convert latitude and longitude to UTM. All coordinates are
    339 #       assumed to be given in the GDA94 datum.
    340 def asc_csiro2sww(bath_dir,
    341                   elevation_dir,
    342                   ucur_dir,
    343                   vcur_dir,
    344                   sww_file,
    345                   minlat=None, maxlat=None,
    346                   minlon=None, maxlon=None,
    347                   zscale=1,
    348                   mean_stage=0,
    349                   fail_on_NaN=True,
    350                   elevation_NaN_filler=0,
    351                   bath_prefix='ba',
    352                   elevation_prefix='el',
    353                   verbose=False):
    354     """
    355     Produce an sww boundary file, from esri ascii data from CSIRO.
    356 
    357     Also convert latitude and longitude to UTM. All coordinates are
    358     assumed to be given in the GDA94 datum.
    359 
    360     assume:
    361     All files are in esri ascii format
    362 
    363     4 types of information
    364     bathymetry
    365     elevation
    366     u velocity
    367     v velocity
    368 
    369     Assumptions
    370     The metadata of all the files is the same
    371     Each type is in a seperate directory
    372     One bath file with extention .000
    373     The time period is less than 24hrs and uniform.
    374     """
    375 
    376     from Scientific.IO.NetCDF import NetCDFFile
    377 
    378     from anuga.coordinate_transforms.redfearn import redfearn
    379 
    380     # So if we want to change the precision it's done here
    381     precision = netcdf_float
    382 
    383     # go in to the bath dir and load the only file,
    384     bath_files = os.listdir(bath_dir)
    385     bath_file = bath_files[0]
    386     bath_dir_file =  bath_dir + os.sep + bath_file
    387     bath_metadata, bath_grid =  _read_asc(bath_dir_file)
    388 
    389     #Use the date.time of the bath file as a basis for
    390     #the start time for other files
    391     base_start = bath_file[-12:]
    392 
    393     #go into the elevation dir and load the 000 file
    394     elevation_dir_file = elevation_dir  + os.sep + elevation_prefix \
    395                          + base_start
    396 
    397     elevation_files = os.listdir(elevation_dir)
    398     ucur_files = os.listdir(ucur_dir)
    399     vcur_files = os.listdir(vcur_dir)
    400     elevation_files.sort()
    401 
    402     # the first elevation file should be the
    403     # file with the same base name as the bath data
    404     assert elevation_files[0] == 'el' + base_start
    405 
    406     number_of_latitudes = bath_grid.shape[0]
    407     number_of_longitudes = bath_grid.shape[1]
    408     number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
    409 
    410     longitudes = [bath_metadata['xllcorner'] + x*bath_metadata['cellsize'] \
    411                   for x in range(number_of_longitudes)]
    412     latitudes = [bath_metadata['yllcorner'] + y*bath_metadata['cellsize'] \
    413                  for y in range(number_of_latitudes)]
    414 
    415      # reverse order of lat, so the first lat represents the first grid row
    416     latitudes.reverse()
    417 
    418     kmin, kmax, lmin, lmax = get_min_max_indices(latitudes[:],longitudes[:],
    419                                                   minlat=minlat, maxlat=maxlat,
    420                                                   minlon=minlon, maxlon=maxlon)
    421 
    422     bath_grid = bath_grid[kmin:kmax,lmin:lmax]
    423     latitudes = latitudes[kmin:kmax]
    424     longitudes = longitudes[lmin:lmax]
    425     number_of_latitudes = len(latitudes)
    426     number_of_longitudes = len(longitudes)
    427     number_of_times = len(os.listdir(elevation_dir))
    428     number_of_points = number_of_latitudes * number_of_longitudes
    429     number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
    430 
    431     #Work out the times
    432     if len(elevation_files) > 1:
    433         # Assume: The time period is less than 24hrs.
    434         time_period = (int(elevation_files[1][-3:]) \
    435                        - int(elevation_files[0][-3:])) * 60*60
    436         times = [x*time_period for x in range(len(elevation_files))]
    437     else:
    438         times = [0.0]
    439 
    440     if verbose:
    441         log.critical('------------------------------------------------')
    442         log.critical('Statistics:')
    443         log.critical('  Extent (lat/lon):')
    444         log.critical('    lat in [%f, %f], len(lat) == %d'
    445                      % (min(latitudes), max(latitudes), len(latitudes)))
    446         log.critical('    lon in [%f, %f], len(lon) == %d'
    447                      % (min(longitudes), max(longitudes), len(longitudes)))
    448         log.critical('    t in [%f, %f], len(t) == %d'
    449                      % (min(times), max(times), len(times)))
    450 
    451     ######### WRITE THE SWW FILE #############
    452 
    453     # NetCDF file definition
    454     outfile = NetCDFFile(sww_file, netcdf_mode_w)
    455 
    456     #Create new file
    457     outfile.institution = 'Geoscience Australia'
    458     outfile.description = 'Converted from XXX'
    459 
    460     #For sww compatibility
    461     outfile.smoothing = 'Yes'
    462     outfile.order = 1
    463 
    464     #Start time in seconds since the epoch (midnight 1/1/1970)
    465     outfile.starttime = starttime = times[0]
    466 
    467     # dimension definitions
    468     outfile.createDimension('number_of_volumes', number_of_volumes)
    469     outfile.createDimension('number_of_vertices', 3)
    470     outfile.createDimension('number_of_points', number_of_points)
    471     outfile.createDimension('number_of_timesteps', number_of_times)
    472 
    473     # variable definitions
    474     outfile.createVariable('x', precision, ('number_of_points',))
    475     outfile.createVariable('y', precision, ('number_of_points',))
    476     outfile.createVariable('elevation', precision, ('number_of_points',))
    477 
    478     #FIXME: Backwards compatibility
    479     #outfile.createVariable('z', precision, ('number_of_points',))
    480     #################################
    481 
    482     outfile.createVariable('volumes', netcdf_int, ('number_of_volumes',
    483                                                    'number_of_vertices'))
    484 
    485     outfile.createVariable('time', precision, ('number_of_timesteps',))
    486 
    487     outfile.createVariable('stage', precision, ('number_of_timesteps',
    488                                                 'number_of_points'))
    489 
    490     outfile.createVariable('xmomentum', precision, ('number_of_timesteps',
    491                                                     'number_of_points'))
    492 
    493     outfile.createVariable('ymomentum', precision, ('number_of_timesteps',
    494                                                     'number_of_points'))
    495 
    496     #Store
    497     from anuga.coordinate_transforms.redfearn import redfearn
    498 
    499     x = num.zeros(number_of_points, num.float)  #Easting
    500     y = num.zeros(number_of_points, num.float)  #Northing
    501 
    502     if verbose: log.critical('Making triangular grid')
    503 
    504     #Get zone of 1st point.
    505     refzone, _, _ = redfearn(latitudes[0], longitudes[0])
    506 
    507     vertices = {}
    508     i = 0
    509     for k, lat in enumerate(latitudes):
    510         for l, lon in enumerate(longitudes):
    511             vertices[l,k] = i
    512 
    513             zone, easting, northing = redfearn(lat, lon)
    514 
    515             #msg = 'Zone boundary crossed at longitude =', lon
    516             #assert zone == refzone, msg
    517             #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
    518             x[i] = easting
    519             y[i] = northing
    520             i += 1
    521 
    522     #Construct 2 triangles per 'rectangular' element
    523     volumes = []
    524     for l in range(number_of_longitudes-1):    #X direction
    525         for k in range(number_of_latitudes-1): #Y direction
    526             v1 = vertices[l,k+1]
    527             v2 = vertices[l,k]
    528             v3 = vertices[l+1,k+1]
    529             v4 = vertices[l+1,k]
    530 
    531             #Note, this is different to the ferrit2sww code
    532             #since the order of the lats is reversed.
    533             volumes.append([v1,v3,v2]) #Upper element
    534             volumes.append([v4,v2,v3]) #Lower element
    535 
    536     volumes = num.array(volumes, num.int)      #array default#
    537 
    538     geo_ref = Geo_reference(refzone, min(x), min(y))
    539     geo_ref.write_NetCDF(outfile)
    540 
    541     # This will put the geo ref in the middle
    542     #geo_ref = Geo_reference(refzone, (max(x)+min(x))/2., (max(x)+min(y))/2.)
    543 
    544     if verbose:
    545         log.critical('------------------------------------------------')
    546         log.critical('More Statistics:')
    547         log.critical('  Extent (/lon):')
    548         log.critical('    x in [%f, %f], len(lat) == %d'
    549                      % (min(x), max(x), len(x)))
    550         log.critical('    y in [%f, %f], len(lon) == %d'
    551                      % (min(y), max(y), len(y)))
    552         log.critical('geo_ref: ', geo_ref)
    553 
    554     z = num.resize(bath_grid,outfile.variables['elevation'][:].shape)
    555     outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
    556     outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
    557     # FIXME (Ole): Remove once viewer has been recompiled and changed
    558     #              to use elevation instead of z
    559     #outfile.variables['z'][:] = z
    560     outfile.variables['elevation'][:] = z
    561     outfile.variables['volumes'][:] = volumes.astype(num.int32) # On Opteron 64
    562 
    563     stage = outfile.variables['stage']
    564     xmomentum = outfile.variables['xmomentum']
    565     ymomentum = outfile.variables['ymomentum']
    566 
    567     outfile.variables['time'][:] = times   #Store time relative
    568 
    569     if verbose: log.critical('Converting quantities')
    570 
    571     n = number_of_times
    572     for j in range(number_of_times):
    573         # load in files
    574         elevation_meta, elevation_grid = \
    575             _read_asc(elevation_dir + os.sep + elevation_files[j])
    576 
    577         _, u_momentum_grid = _read_asc(ucur_dir + os.sep + ucur_files[j])
    578         _, v_momentum_grid = _read_asc(vcur_dir + os.sep + vcur_files[j])
    579 
    580         #cut matrix to desired size
    581         elevation_grid = elevation_grid[kmin:kmax,lmin:lmax]
    582         u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
    583         v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
    584 
    585         # handle missing values
    586         missing = (elevation_grid == elevation_meta['NODATA_value'])
    587         if num.sometrue (missing):
    588             if fail_on_NaN:
    589                 msg = 'File %s contains missing values' \
    590                       % (elevation_files[j])
    591                 raise DataMissingValuesError, msg
    592             else:
    593                 elevation_grid = elevation_grid*(missing==0) \
    594                                  + missing*elevation_NaN_filler
    595 
    596         if verbose and j % ((n+10)/10) == 0: log.critical('  Doing %d of %d'
    597                                                           % (j, n))
    598 
    599         i = 0
    600         for k in range(number_of_latitudes):      #Y direction
    601             for l in range(number_of_longitudes): #X direction
    602                 w = zscale*elevation_grid[k,l] + mean_stage
    603                 stage[j,i] = w
    604                 h = w - z[i]
    605                 xmomentum[j,i] = u_momentum_grid[k,l]*h
    606                 ymomentum[j,i] = v_momentum_grid[k,l]*h
    607                 i += 1
    608 
    609     outfile.close()
    610 
    611 
    612 
    613 def _read_asc(filename, verbose=False):
    614     """Read esri file from the following ASCII format (.asc)
    615 
    616     Example:
    617 
    618     ncols         3121
    619     nrows         1800
    620     xllcorner     722000
    621     yllcorner     5893000
    622     cellsize      25
    623     NODATA_value  -9999
    624     138.3698 137.4194 136.5062 135.5558 ..........
    625 
    626     filename Path to the file to read.
    627     verbose True if this function is to be verbose.
    628     """
    629 
    630     datafile = open(filename)
    631 
    632     if verbose: log.critical('Reading DEM from %s' % filename)
    633 
    634     lines = datafile.readlines()
    635     datafile.close()
    636 
    637     if verbose: log.critical('Got %d lines' % len(lines))
    638 
    639     ncols = int(lines.pop(0).split()[1].strip())
    640     nrows = int(lines.pop(0).split()[1].strip())
    641     xllcorner = float(lines.pop(0).split()[1].strip())
    642     yllcorner = float(lines.pop(0).split()[1].strip())
    643     cellsize = float(lines.pop(0).split()[1].strip())
    644     NODATA_value = float(lines.pop(0).split()[1].strip())
    645 
    646     assert len(lines) == nrows
    647 
    648     #Store data
    649     grid = []
    650 
    651     n = len(lines)
    652     for i, line in enumerate(lines):
    653         cells = line.split()
    654         assert len(cells) == ncols
    655         grid.append(num.array([float(x) for x in cells]))
    656     grid = num.array(grid)
    657 
    658     return {'xllcorner':xllcorner,
    659             'yllcorner':yllcorner,
    660             'cellsize':cellsize,
    661             'NODATA_value':NODATA_value}, grid
    662318
    663319
  • trunk/anuga_core/source/anuga/file_conversion/test_file_conversion.py

    r7744 r7758  
    1111from anuga.config import netcdf_float, epsilon, g
    1212from Scientific.IO.NetCDF import NetCDFFile
     13from anuga.file_conversion.file_conversion import tsh2sww, \
     14                        pmesh_to_domain_instance
     15
    1316import sys
    1417import unittest
Note: See TracChangeset for help on using the changeset viewer.