Changeset 7767


Ignore:
Timestamp:
Jun 2, 2010, 3:00:14 PM (14 years ago)
Author:
hudson
Message:

Added sww module.

Location:
trunk/anuga_core/source/anuga/file
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/file/sww.py

    r7762 r7767  
    984984        log.critical('------------------------------------------------')
    985985
     986
     987
     988
     989##
     990# @brief Get the extents of a NetCDF data file.
     991# @param file_name The path to the SWW file.
     992# @return A list of x, y, z and stage limits (min, max).
     993def extent_sww(file_name):
     994    """Read in an sww file.
     995
     996    Input:
     997    file_name - the sww file
     998
     999    Output:
     1000    A list: [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
     1001    """
     1002
     1003    from Scientific.IO.NetCDF import NetCDFFile
     1004
     1005    #Get NetCDF
     1006    fid = NetCDFFile(file_name, netcdf_mode_r)
     1007
     1008    # Get the variables
     1009    x = fid.variables['x'][:]
     1010    y = fid.variables['y'][:]
     1011    stage = fid.variables['stage'][:]
     1012
     1013    fid.close()
     1014
     1015    return [min(x), max(x), min(y), max(y), num.min(stage), num.max(stage)]
     1016
     1017
     1018##
     1019# @brief
     1020# @param filename
     1021# @param boundary
     1022# @param t
     1023# @param fail_if_NaN
     1024# @param NaN_filler
     1025# @param verbose
     1026# @param very_verbose
     1027# @return
     1028def load_sww_as_domain(filename, boundary=None, t=None,
     1029               fail_if_NaN=True, NaN_filler=0,
     1030               verbose=False, very_verbose=False):
     1031    """
     1032    Usage: domain = load_sww_as_domain('file.sww',
     1033                        t=time (default = last time in file))
     1034
     1035    Boundary is not recommended if domain.smooth is not selected, as it
     1036    uses unique coordinates, but not unique boundaries. This means that
     1037    the boundary file will not be compatable with the coordinates, and will
     1038    give a different final boundary, or crash.
     1039    """
     1040   
     1041    from Scientific.IO.NetCDF import NetCDFFile
     1042    from shallow_water import Domain
     1043
     1044    # initialise NaN.
     1045    NaN = 9.969209968386869e+036
     1046
     1047    if verbose: log.critical('Reading from %s' % filename)
     1048
     1049    fid = NetCDFFile(filename, netcdf_mode_r)    # Open existing file for read
     1050    time = fid.variables['time']       # Timesteps
     1051    if t is None:
     1052        t = time[-1]
     1053    time_interp = get_time_interp(time,t)
     1054
     1055    # Get the variables as numeric arrays
     1056    x = fid.variables['x'][:]                   # x-coordinates of vertices
     1057    y = fid.variables['y'][:]                   # y-coordinates of vertices
     1058    elevation = fid.variables['elevation']      # Elevation
     1059    stage = fid.variables['stage']              # Water level
     1060    xmomentum = fid.variables['xmomentum']      # Momentum in the x-direction
     1061    ymomentum = fid.variables['ymomentum']      # Momentum in the y-direction
     1062
     1063    starttime = fid.starttime[0]
     1064    volumes = fid.variables['volumes'][:]       # Connectivity
     1065    coordinates = num.transpose(num.asarray([x.tolist(), y.tolist()]))
     1066    # FIXME (Ole): Something like this might be better:
     1067    #                 concatenate((x, y), axis=1)
     1068    # or              concatenate((x[:,num.newaxis], x[:,num.newaxis]), axis=1)
     1069
     1070    conserved_quantities = []
     1071    interpolated_quantities = {}
     1072    other_quantities = []
     1073
     1074    # get geo_reference
     1075    try:                             # sww files don't have to have a geo_ref
     1076        geo_reference = Geo_reference(NetCDFObject=fid)
     1077    except: # AttributeError, e:
     1078        geo_reference = None
     1079
     1080    if verbose: log.critical('    getting quantities')
     1081
     1082    for quantity in fid.variables.keys():
     1083        dimensions = fid.variables[quantity].dimensions
     1084        if 'number_of_timesteps' in dimensions:
     1085            conserved_quantities.append(quantity)
     1086            interpolated_quantities[quantity] = \
     1087                  interpolated_quantity(fid.variables[quantity][:], time_interp)
     1088        else:
     1089            other_quantities.append(quantity)
     1090
     1091    other_quantities.remove('x')
     1092    other_quantities.remove('y')
     1093    #other_quantities.remove('z')
     1094    other_quantities.remove('volumes')
     1095    try:
     1096        other_quantities.remove('stage_range')
     1097        other_quantities.remove('xmomentum_range')
     1098        other_quantities.remove('ymomentum_range')
     1099        other_quantities.remove('elevation_range')
     1100    except:
     1101        pass
     1102
     1103    conserved_quantities.remove('time')
     1104
     1105    if verbose: log.critical('    building domain')
     1106
     1107    #    From domain.Domain:
     1108    #    domain = Domain(coordinates, volumes,\
     1109    #                    conserved_quantities = conserved_quantities,\
     1110    #                    other_quantities = other_quantities,zone=zone,\
     1111    #                    xllcorner=xllcorner, yllcorner=yllcorner)
     1112
     1113    # From shallow_water.Domain:
     1114    coordinates = coordinates.tolist()
     1115    volumes = volumes.tolist()
     1116    # FIXME:should this be in mesh? (peter row)
     1117    if fid.smoothing == 'Yes':
     1118        unique = False
     1119    else:
     1120        unique = True
     1121    if unique:
     1122        coordinates, volumes, boundary = weed(coordinates, volumes,boundary)
     1123
     1124     
     1125   
     1126    try:
     1127        domain = Domain(coordinates, volumes, boundary)
     1128    except AssertionError, e:
     1129        fid.close()
     1130        msg = 'Domain could not be created: %s. ' \
     1131              'Perhaps use "fail_if_NaN=False and NaN_filler = ..."' % e
     1132        raise DataDomainError, msg
     1133
     1134    if not boundary is None:
     1135        domain.boundary = boundary
     1136
     1137    domain.geo_reference = geo_reference
     1138
     1139    domain.starttime = float(starttime) + float(t)
     1140    domain.time = 0.0
     1141
     1142    for quantity in other_quantities:
     1143        try:
     1144            NaN = fid.variables[quantity].missing_value
     1145        except:
     1146            pass                       # quantity has no missing_value number
     1147        X = fid.variables[quantity][:]
     1148        if very_verbose:
     1149            log.critical('       %s' % str(quantity))
     1150            log.critical('        NaN = %s' % str(NaN))
     1151            log.critical('        max(X)')
     1152            log.critical('       %s' % str(max(X)))
     1153            log.critical('        max(X)==NaN')
     1154            log.critical('       %s' % str(max(X)==NaN))
     1155            log.critical('')
     1156        if max(X) == NaN or min(X) == NaN:
     1157            if fail_if_NaN:
     1158                msg = 'quantity "%s" contains no_data entry' % quantity
     1159                raise DataMissingValuesError, msg
     1160            else:
     1161                data = (X != NaN)
     1162                X = (X*data) + (data==0)*NaN_filler
     1163        if unique:
     1164            X = num.resize(X, (len(X)/3, 3))
     1165        domain.set_quantity(quantity, X)
     1166    #
     1167    for quantity in conserved_quantities:
     1168        try:
     1169            NaN = fid.variables[quantity].missing_value
     1170        except:
     1171            pass                       # quantity has no missing_value number
     1172        X = interpolated_quantities[quantity]
     1173        if very_verbose:
     1174            log.critical('       %s' % str(quantity))
     1175            log.critical('        NaN = %s' % str(NaN))
     1176            log.critical('        max(X)')
     1177            log.critical('       %s' % str(max(X)))
     1178            log.critical('        max(X)==NaN')
     1179            log.critical('       %s' % str(max(X)==NaN))
     1180            log.critical('')
     1181        if max(X) == NaN or min(X) == NaN:
     1182            if fail_if_NaN:
     1183                msg = 'quantity "%s" contains no_data entry' % quantity
     1184                raise DataMissingValuesError, msg
     1185            else:
     1186                data = (X != NaN)
     1187                X = (X*data) + (data==0)*NaN_filler
     1188        if unique:
     1189            X = num.resize(X, (X.shape[0]/3, 3))
     1190        domain.set_quantity(quantity, X)
     1191
     1192    fid.close()
     1193
     1194    return domain
     1195
     1196
     1197##
     1198# @brief
     1199# @parm time
     1200# @param t
     1201# @return An (index, ration) tuple.
     1202def get_time_interp(time, t=None):
     1203    #Finds the ratio and index for time interpolation.
     1204    #It is borrowed from previous abstract_2d_finite_volumes code.
     1205    if t is None:
     1206        t=time[-1]
     1207        index = -1
     1208        ratio = 0.
     1209    else:
     1210        T = time
     1211        tau = t
     1212        index=0
     1213        msg = 'Time interval derived from file %s [%s:%s]' \
     1214              % ('FIXMEfilename', T[0], T[-1])
     1215        msg += ' does not match model time: %s' % tau
     1216        if tau < time[0]: raise DataTimeError, msg
     1217        if tau > time[-1]: raise DataTimeError, msg
     1218        while tau > time[index]: index += 1
     1219        while tau < time[index]: index -= 1
     1220        if tau == time[index]:
     1221            #Protect against case where tau == time[-1] (last time)
     1222            # - also works in general when tau == time[i]
     1223            ratio = 0
     1224        else:
     1225            #t is now between index and index+1
     1226            ratio = (tau - time[index])/(time[index+1] - time[index])
     1227
     1228    return (index, ratio)
     1229
     1230
     1231
     1232##
     1233# @brief Interpolate a quantity wrt time.
     1234# @param saved_quantity The quantity to interpolate.
     1235# @param time_interp (index, ratio)
     1236# @return A vector of interpolated values.
     1237def interpolated_quantity(saved_quantity, time_interp):
     1238    '''Given an index and ratio, interpolate quantity with respect to time.'''
     1239
     1240    index, ratio = time_interp
     1241
     1242    Q = saved_quantity
     1243
     1244    if ratio > 0:
     1245        q = (1-ratio)*Q[index] + ratio*Q[index+1]
     1246    else:
     1247        q = Q[index]
     1248
     1249    #Return vector of interpolated values
     1250    return q
     1251
     1252
     1253
     1254
     1255##
     1256# @brief
     1257# @param coordinates
     1258# @param volumes
     1259# @param boundary
     1260def weed(coordinates, volumes, boundary=None):
     1261    if isinstance(coordinates, num.ndarray):
     1262        coordinates = coordinates.tolist()
     1263    if isinstance(volumes, num.ndarray):
     1264        volumes = volumes.tolist()
     1265
     1266    unique = False
     1267    point_dict = {}
     1268    same_point = {}
     1269    for i in range(len(coordinates)):
     1270        point = tuple(coordinates[i])
     1271        if point_dict.has_key(point):
     1272            unique = True
     1273            same_point[i] = point
     1274            #to change all point i references to point j
     1275        else:
     1276            point_dict[point] = i
     1277            same_point[i] = point
     1278
     1279    coordinates = []
     1280    i = 0
     1281    for point in point_dict.keys():
     1282        point = tuple(point)
     1283        coordinates.append(list(point))
     1284        point_dict[point] = i
     1285        i += 1
     1286
     1287    for volume in volumes:
     1288        for i in range(len(volume)):
     1289            index = volume[i]
     1290            if index > -1:
     1291                volume[i] = point_dict[same_point[index]]
     1292
     1293    new_boundary = {}
     1294    if not boundary is None:
     1295        for segment in boundary.keys():
     1296            point0 = point_dict[same_point[segment[0]]]
     1297            point1 = point_dict[same_point[segment[1]]]
     1298            label = boundary[segment]
     1299            #FIXME should the bounday attributes be concaterated
     1300            #('exterior, pond') or replaced ('pond')(peter row)
     1301
     1302            if new_boundary.has_key((point0, point1)):
     1303                new_boundary[(point0,point1)] = new_boundary[(point0, point1)]
     1304
     1305            elif new_boundary.has_key((point1, point0)):
     1306                new_boundary[(point1,point0)] = new_boundary[(point1, point0)]
     1307            else: new_boundary[(point0, point1)] = label
     1308
     1309        boundary = new_boundary
     1310
     1311    return coordinates, volumes, boundary
Note: See TracChangeset for help on using the changeset viewer.