Changeset 7688


Ignore:
Timestamp:
Apr 19, 2010, 3:11:12 PM (14 years ago)
Author:
fountain
Message:

updates to bunbury storm surge model

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py

    r7684 r7688  
    4242print "Elevation file name: ", elevation_file_name
    4343
    44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780403.18*.gz'))
     44grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780404.0*.gz'))
    4545grid_file_names.sort()
    46 number_of_timesteps = len(grid_file_names)
     46
     47# grid_file_names_temp = []
     48# for i in range(0, 120, 4):
     49    # grid_file_names_temp.append(grid_file_names[i])
     50
     51# grid_file_names = grid_file_names_temp
     52   
     53# number_of_timesteps = len(grid_file_names)
    4754
    4855print "Number of timesteps: ", number_of_timesteps
     
    139146number_of_points = ncols * nrows
    140147
    141 # insert nodata test here
     148# ---------------------------------------------
     149# Create mask of no_data values across stage, speed and theta to ensure all three quantities
     150# have no_data values in corresponding cells - this is a work-around until GEMS can produce a
     151# a consistant dataset
     152# ---------------------------------------------
    142153
    143154for i in xrange(number_of_timesteps):
    144155    j = number_of_points * i
    145156    k = j + number_of_points
    146     d = stage[j:k] - elevation     # Assumes elevations below sea level are negative
    147     depth[j:k] = [0 if x<-9000 else x for x in d]   # sets stage to zero for no-data values of stage
     157    # no_value_mask = (stage[j:k] < -9000) + (speed[j:k] < -9000) + (theta[j:k] < -9000)
     158    no_value_index = numpy.where(((stage[j:k] < -100) + (speed[j:k] < 0) + (theta[j:k] < 0)) == True)[0]
     159    depth[j:k] = stage[j:k] - elevation
     160    numpy.put(stage[j:k], no_value_index, -9999)
     161    numpy.put(speed[j:k], no_value_index, -9999)
     162    numpy.put(theta[j:k], no_value_index, -9999)
     163    numpy.put(depth[j:k], no_value_index, 0)
     164
     165# for i in xrange(number_of_timesteps):
     166    # j = number_of_points * i
     167    # k = j + number_of_points
     168    # d = stage[j:k] - elevation     # Assumes elevations below sea level are negative
     169    # depth[j:k] = [0 if x<-9000 else x for x in d]   # sets stage to zero for no-data values of stage
    148170   
    149171    # depth[j:k] = [x if x>=0 and elevation>0 else 0 for x in d] # This is to correct for GEMS data which assumes stage=0 onshore when dry but should = elevation
     
    160182
    161183for i, t in enumerate(theta):
     184   
     185    mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north
     186    my = momentum[i]*math.cos(math.radians(t))
    162187
    163188    if t > 0.0 and t < 90.0:
    164         mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north
    165         my = momentum[i]*math.cos(math.radians(t))
    166189        assert mx > 0 and my > 0
     190       
    167191    elif t > 90.0 and t < 180.0:
    168         mx = momentum[i]*math.sin(math.radians(t))
    169         my = momentum[i]*math.cos(math.radians(t))
    170192        assert mx > 0 and my < 0
     193       
    171194    elif t > 180.0 and t < 270.0:
    172         mx = momentum[i]*math.cos(math.radians(t))
    173         my = momentum[i]*math.sin(math.radians(t))
    174195        assert mx < 0 and my < 0
     196       
    175197    elif t > 270.0 and t < 360.0:
    176         mx = momentum[i]*math.sin(math.radians(t))
    177         my = momentum[i]*math.cos(math.radians(t))
    178198        assert mx < 0 and my > 0
     199       
    179200    elif t == 0.0 or t == 360.0:
    180         mx = 0
    181         my = momentum[i]
    182         assert my >= 0
     201        assert my > 0
     202       
    183203    elif t == 90.0:
    184         mx = momentum[i]
    185         my = 0
    186         assert mx >= 0
     204        assert mx > 0
     205       
    187206    elif t == 180.0:
    188         mx = 0
    189         my = -momentum[i]
    190         assert my <= 0
     207        assert my < 0
     208       
    191209    elif t == 270.0:
    192         mx = -momentum[i]
    193         my = 0
    194         assert mx <= 0
     210        assert mx < 0
     211       
    195212    elif t == -9999:
    196213        mx = 0
     
    219236assert momentum.size == stage.size
    220237
    221 #-----------------------------
     238# -----------------------------
    222239# Create the STS file
    223 #-----------------------------
     240# -----------------------------
    224241print "Creating the STS NetCDF file"
    225242
Note: See TracChangeset for help on using the changeset viewer.