Changeset 7688
- Timestamp:
- Apr 19, 2010, 3:11:12 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py
r7684 r7688 42 42 print "Elevation file name: ", elevation_file_name 43 43 44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt1978040 3.18*.gz'))44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780404.0*.gz')) 45 45 grid_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) 47 54 48 55 print "Number of timesteps: ", number_of_timesteps … … 139 146 number_of_points = ncols * nrows 140 147 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 # --------------------------------------------- 142 153 143 154 for i in xrange(number_of_timesteps): 144 155 j = number_of_points * i 145 156 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 148 170 149 171 # 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 … … 160 182 161 183 for 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)) 162 187 163 188 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 north165 my = momentum[i]*math.cos(math.radians(t))166 189 assert mx > 0 and my > 0 190 167 191 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))170 192 assert mx > 0 and my < 0 193 171 194 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))174 195 assert mx < 0 and my < 0 196 175 197 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))178 198 assert mx < 0 and my > 0 199 179 200 elif t == 0.0 or t == 360.0: 180 mx = 0 181 my = momentum[i] 182 assert my >= 0 201 assert my > 0 202 183 203 elif t == 90.0: 184 mx = momentum[i] 185 my = 0 186 assert mx >= 0 204 assert mx > 0 205 187 206 elif t == 180.0: 188 mx = 0 189 my = -momentum[i] 190 assert my <= 0 207 assert my < 0 208 191 209 elif t == 270.0: 192 mx = -momentum[i] 193 my = 0 194 assert mx <= 0 210 assert mx < 0 211 195 212 elif t == -9999: 196 213 mx = 0 … … 219 236 assert momentum.size == stage.size 220 237 221 # -----------------------------238 # ----------------------------- 222 239 # Create the STS file 223 # -----------------------------240 # ----------------------------- 224 241 print "Creating the STS NetCDF file" 225 242
Note: See TracChangeset
for help on using the changeset viewer.