Changeset 3964


Ignore:
Timestamp:
Nov 10, 2006, 3:47:17 PM (17 years ago)
Author:
duncan
Message:

comments and minor changes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r3961 r3964  
    41084108    Lat and lon are assumed to be in decimal degrees.
    41094109    NOTE: min Lon is the most east boundary.
     4110   
     4111    origin is a 3-tuple with geo referenced
     4112    UTM coordinates (zone, easting, northing)
     4113    It will be the origin of the sww file. This shouldn't be used,
     4114    since all of anuga should be able to handle an arbitary origin.
     4115
    41104116
    41114117    URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
     
    41384144   
    41394145def urs2nc(basename_in = 'o', basename_out = 'urs'):
     4146    """
     4147    Convert the 3 urs files to 4 nc files.
     4148
     4149    The name of the urs file names must be;
     4150    [basename_in]-z-mux
     4151    [basename_in]-e-mux
     4152    [basename_in]-n-mux
     4153   
     4154    """
    41404155    files_in = [basename_in+'-z-mux',
    41414156                basename_in+'-e-mux',
     
    41584173                                         file_out,
    41594174                                         quantity)
    4160         #print "lon",lon
    4161         #print "lat",lat
    41624175        if hashed_elevation == None:
    41634176            elevation_file = basename_out+'_e.nc'
    4164             write_elevation_sww(elevation_file,
     4177            write_elevation_nc(elevation_file,
    41654178                                lon,
    41664179                                lat,
     
    41754188def _binary_c2nc(file_in, file_out, quantity):
    41764189    """
    4177 
    4178     return the depth info, so it can be written to a file
     4190    Reads in a quantity urs file and writes a quantity nc file.
     4191    additionally, returns the depth and lat, long info,
     4192    so it can be written to a file.
    41794193    """
    41804194    columns = 3 # long, lat , depth
    4181     #FIXME use mux_file, not f
    4182     f = open(file_in, 'rb')
     4195    mux_file = open(file_in, 'rb')
    41834196   
    41844197    # Number of points/stations
    4185     (points_num,)= unpack('i',f.read(4))
    4186 
    4187     #print 'points_num', points_num
    4188     #import sys; sys.exit()
     4198    (points_num,)= unpack('i',mux_file.read(4))
     4199
    41894200    # nt, int - Number of time steps
    4190     (time_step_count,)= unpack('i',f.read(4))
     4201    (time_step_count,)= unpack('i',mux_file.read(4))
    41914202
    41924203    #dt, float - time step, seconds
    4193     (time_step,) = unpack('f', f.read(4))
     4204    (time_step,) = unpack('f', mux_file.read(4))
    41944205   
    41954206    msg = "Bad data in the mux file."
    41964207    if points_num < 0:
    4197         f.close()
     4208        mux_file.close()
    41984209        raise ANUGAError, msg
    41994210    if time_step_count < 0:
    4200         f.close()
     4211        mux_file.close()
    42014212        raise ANUGAError, msg
    42024213    if time_step < 0:
    4203         f.close()
     4214        mux_file.close()
    42044215        raise ANUGAError, msg
    42054216   
    42064217    lonlatdep = p_array.array('f')
    4207     ####points_num = 2914
    4208     lonlatdep.read(f, columns * points_num)
     4218    lonlatdep.read(mux_file, columns * points_num)
    42094219    lonlatdep = array(lonlatdep, typecode=Float)   
    42104220    lonlatdep = reshape(lonlatdep, (points_num, columns))
    4211    
    4212     #print "lonlatdep", lonlatdep
    42134221   
    42144222    lon, lat, depth = lon_lat2grid(lonlatdep)
     
    42164224    lon_sorted.sort()
    42174225
    4218     if not allclose(lon, lon_sorted):
     4226    if not lon == lon_sorted:
    42194227        msg = "Longitudes in mux file are not in ascending order"
    42204228        raise IOError, msg
     
    42224230    lat_sorted.sort()
    42234231
    4224     if not allclose(lat, lat_sorted):
     4232    if not lat == lat_sorted:
    42254233        msg = "Latitudes in mux file are not in ascending order"
    42264234   
     
    42354243        #Read in a time slice  from mux file 
    42364244        hz_p_array = p_array.array('f')
    4237         hz_p_array.read(f, points_num)
     4245        hz_p_array.read(mux_file, points_num)
    42384246        hz_p = array(hz_p_array, typecode=Float)
    42394247        hz_p = reshape(hz_p, (len(lon), len(lat)))
     
    42424250        #write time slice to nc file
    42434251        nc_file.store_timestep(hz_p)
    4244     #FIXME should I close the mux file here?
    4245     f.close()
     4252    mux_file.close()
    42464253    nc_file.close()
    42474254
     
    42494256   
    42504257
    4251 def write_elevation_sww(file_out, lon, lat, depth_vector):
    4252      
     4258def write_elevation_nc(file_out, lon, lat, depth_vector):
     4259    """
     4260    Write an nc elevation file.
     4261    """
     4262   
    42534263    # NetCDF file definition
    42544264    outfile = NetCDFFile(file_out, 'w')
     
    42674277
    42684278    depth = reshape(depth_vector, ( len(lat), len(lon)))
    4269     #print "depth",depth
    42704279    outfile.variables[zname][:]= depth
    42714280   
     
    42734282   
    42744283def nc_lon_lat_header(outfile, lon, lat):
     4284    """
     4285    outfile is the netcdf file handle.
     4286    lon - a list/array of the longitudes
     4287    lat - a list/array of the latitudes
     4288    """
    42754289   
    42764290    outfile.institution = 'Geoscience Australia'
     
    43124326    num_points = long_lat_dep.shape[0]
    43134327    this_rows_long = long_lat_dep[0,LONG]
    4314    
     4328
    43154329    # Count the length of unique latitudes
    43164330    i = 0
    43174331    while long_lat_dep[i,LONG] == this_rows_long and i < num_points:
    43184332        i += 1
    4319          
     4333    # determine the lats and longsfrom the grid
    43204334    lat = long_lat_dep[:i, LAT]       
    4321     long = long_lat_dep[::i, LONG]     
     4335    long = long_lat_dep[::i, LONG]
     4336   
    43224337    lenlong = len(long)
    43234338    lenlat = len(lat)
    4324 #    print 'len lat', lat, len(lat)
    4325 #    print 'len long', long, len(long)
     4339    #print 'len lat', lat, len(lat)
     4340    #print 'len long', long, len(long)
    43264341         
    43274342    msg = 'Input data is not gridded'     
Note: See TracChangeset for help on using the changeset viewer.