Changeset 7766 for trunk/anuga_core/source/anuga/file
- Timestamp:
- Jun 2, 2010, 1:15:38 PM (15 years ago)
- Location:
- trunk/anuga_core/source/anuga/file
- Files:
-
- 1 added
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/file/mux.py
r7760 r7766 1 """ 2 Read a mux2 file. 3 """ 4 1 5 from anuga.utilities.numerical_tools import ensure_numeric 2 6 import numpy as num 3 7 4 8 ################################################################################ 5 9 # READ MUX2 FILES line of points 6 10 ################################################################################ 11 12 WAVEHEIGHT_MUX_LABEL = '-z-mux' 13 EAST_VELOCITY_LABEL = '-e-mux' 14 NORTH_VELOCITY_LABEL = '-n-mux' 7 15 8 16 WAVEHEIGHT_MUX2_LABEL = '-z-mux2' … … 43 51 # Convert verbose to int C flag 44 52 if verbose: 45 verbose =153 verbose = 1 46 54 else: 47 verbose =055 verbose = 0 48 56 49 57 if weights is None: … … 102 110 longitudes = num.zeros(number_of_selected_stations, num.float) 103 111 elevation = num.zeros(number_of_selected_stations, num.float) 104 quantity = num.zeros((number_of_selected_stations, parameters_index), num.float) 112 quantity = num.zeros((number_of_selected_stations, parameters_index), \ 113 num.float) 105 114 106 115 starttime = 1e16 … … 116 125 117 126 118 ##119 # @brief ??120 # @param mux_times ??121 # @param mint ??122 # @param maxt ??123 # @return ??124 def read_time_from_mux(mux_times, mint, maxt):125 """126 """127 127 128 if mint == None:129 mux_times_start_i = 0130 else:131 mux_times_start_i = num.searchsorted(mux_times, mint)132 133 if maxt == None:134 mux_times_fin_i = len(mux_times)135 else:136 maxt += 0.5 # so if you specify a time where there is137 # data that time will be included138 mux_times_fin_i = num.searchsorted(mux_times, maxt)139 140 return mux_times_start_i, mux_times_fin_i141 142 143 -
trunk/anuga_core/source/anuga/file/sts.py
r7762 r7766 284 284 except: 285 285 msg = 'Cannot open %s' % stsname + '.sts' 286 raise msg286 raise IOError(msg) 287 287 288 288 xllcorner = fid.xllcorner[0] -
trunk/anuga_core/source/anuga/file/test_mux.py
r7762 r7766 10 10 from anuga.coordinate_transforms.geo_reference import Geo_reference 11 11 12 from mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \ 13 NORTH_VELOCITY_LABEL 14 12 15 from mux import WAVEHEIGHT_MUX2_LABEL, EAST_VELOCITY_MUX2_LABEL, \ 13 16 NORTH_VELOCITY_MUX2_LABEL … … 15 18 from mux import read_mux2_py 16 19 from anuga.file_conversion.urs2sts import urs2sts 17 18 class TestCase(unittest.TestCase): 20 from anuga.file.urs import Read_urs 21 22 class Test_Mux(unittest.TestCase): 19 23 def setUp(self): 20 24 pass … … 1487 1491 1488 1492 self.delete_mux(files) 1489 1493 1494 def test_Urs_points(self): 1495 time_step_count = 3 1496 time_step = 2 1497 lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,115)] 1498 base_name, files = self.write_mux(lat_long_points, 1499 time_step_count, time_step) 1500 for file in files: 1501 urs = Read_urs(file) 1502 assert time_step_count == urs.time_step_count 1503 assert time_step == urs.time_step 1504 1505 for lat_lon, dep in map(None, lat_long_points, urs.lonlatdep): 1506 _ , e, n = redfearn(lat_lon[0], lat_lon[1]) 1507 assert num.allclose(n, dep[2]) 1508 1509 count = 0 1510 for slice in urs: 1511 count += 1 1512 #print slice 1513 for lat_lon, quantity in map(None, lat_long_points, slice): 1514 _ , e, n = redfearn(lat_lon[0], lat_lon[1]) 1515 #print "quantity", quantity 1516 #print "e", e 1517 #print "n", n 1518 if file[-5:] == WAVEHEIGHT_MUX_LABEL[-5:] or \ 1519 file[-5:] == NORTH_VELOCITY_LABEL[-5:] : 1520 assert num.allclose(e, quantity) 1521 if file[-5:] == EAST_VELOCITY_LABEL[-5:]: 1522 assert num.allclose(n, quantity) 1523 assert count == time_step_count 1524 1525 self.delete_mux(files) 1526 1527 1528 1490 1529 1491 1530 ################################################################################ 1492 1531 1493 1532 if __name__ == "__main__": 1494 suite = unittest.makeSuite(Test Case,'test')1533 suite = unittest.makeSuite(Test_Mux,'test') 1495 1534 runner = unittest.TextTestRunner() #verbosity=2) 1496 1535 runner.run(suite) -
trunk/anuga_core/source/anuga/file/urs.py
r7760 r7766 1 ## 2 # @brief 1 from struct import pack, unpack 2 import array as p_array 3 import numpy as num 4 5 6 from anuga.coordinate_transforms.geo_reference import Geo_reference 7 8 from anuga.geospatial_data.geospatial_data import ensure_absolute, \ 9 Geospatial_data 10 11 from anuga.coordinate_transforms.redfearn import redfearn 12 3 13 class Read_urs: 4 14 """ … … 88 98 self.mux_file.close() 89 99 100 101 102 103 ### PRODUCING THE POINTS NEEDED FILE ### 104 105 # Ones used for FESA 2007 results 106 #LL_LAT = -50.0 107 #LL_LONG = 80.0 108 #GRID_SPACING = 1.0/60.0 109 #LAT_AMOUNT = 4800 110 #LONG_AMOUNT = 3600 111 112 113 ## 114 # @brief 115 # @param file_name 116 # @param boundary_polygon 117 # @param zone 118 # @param ll_lat 119 # @param ll_long 120 # @param grid_spacing 121 # @param lat_amount 122 # @param long_amount 123 # @param isSouthernHemisphere 124 # @param export_csv 125 # @param use_cache 126 # @param verbose True if this function is to be verbose. 127 # @return 128 129 def save_boundary_as_urs(file_name, boundary_polygon, zone, 130 ll_lat, ll_long, 131 grid_spacing, 132 lat_amount, long_amount, 133 isSouthernHemisphere=True, 134 export_csv=False, use_cache=False, 135 verbose=False): 136 """ 137 Given the info to replicate the URS grid and a polygon output 138 a file that specifies the cloud of boundary points for URS. 139 140 This creates a .urs file. This is in the format used by URS; 141 1st line is the number of points, 142 each line after represents a point,in lats and longs. 143 144 Note: The polygon cannot cross zones or hemispheres. 145 146 A work-a-round for different zones or hemispheres is to run this twice, 147 once for each zone, and then combine the output. 148 149 file_name - name of the urs file produced for David. 150 boundary_polygon - a list of points that describes a polygon. 151 The last point is assumed ot join the first point. 152 This is in UTM (lat long would be better though) 153 154 This is info about the URS model that needs to be inputted. 155 156 ll_lat - lower left latitude, in decimal degrees 157 ll-long - lower left longitude, in decimal degrees 158 grid_spacing - in deciamal degrees 159 lat_amount - number of latitudes 160 long_amount- number of longs 161 162 Don't add the file extension. It will be added. 163 """ 164 165 geo = calculate_boundary_points(boundary_polygon, zone, ll_lat, ll_long, 166 grid_spacing, 167 lat_amount, long_amount, isSouthernHemisphere, 168 use_cache, verbose) 169 170 if not file_name[-4:] == ".urs": 171 file_name += ".urs" 172 173 geo.export_points_file(file_name, isSouthHemisphere=isSouthernHemisphere) 174 175 if export_csv: 176 if file_name[-4:] == ".urs": 177 file_name = file_name[:-4] + ".csv" 178 geo.export_points_file(file_name) 179 180 return geo 181 182 183 ## 184 # @brief 185 # @param boundary_polygon 186 # @param zone 187 # @param ll_lat 188 # @param ll_long 189 # @param grid_spacing 190 # @param lat_amount 191 # @param long_amount 192 # @param isSouthHemisphere 193 # @param use_cache 194 # @param verbose 195 def calculate_boundary_points(boundary_polygon, zone, ll_lat, 196 ll_long, grid_spacing, 197 lat_amount, long_amount, isSouthHemisphere=True, 198 use_cache=False, verbose=False): 199 args = (boundary_polygon, 200 zone, ll_lat, 201 ll_long, grid_spacing, 202 lat_amount, long_amount, isSouthHemisphere) 203 kwargs = {} 204 205 if use_cache is True: 206 try: 207 from anuga.caching import cache 208 except: 209 msg = 'Caching was requested, but caching module' \ 210 'could not be imported' 211 raise msg 212 213 geo = cache(_calculate_boundary_points, 214 args, kwargs, 215 verbose=verbose, 216 compression=False) 217 else: 218 geo = apply(_calculate_boundary_points, args, kwargs) 219 220 return geo 221 222 223 ## 224 # @brief 225 # @param boundary_polygon An iterable of points that describe a polygon. 226 # @param zone 227 # @param ll_lat Lower left latitude, in decimal degrees 228 # @param ll_long Lower left longitude, in decimal degrees 229 # @param grid_spacing Grid spacing in decimal degrees. 230 # @param lat_amount 231 # @param long_amount 232 # @param isSouthHemisphere 233 def _calculate_boundary_points(boundary_polygon, 234 zone, ll_lat, 235 ll_long, grid_spacing, 236 lat_amount, long_amount, 237 isSouthHemisphere): 238 """ 239 boundary_polygon - a list of points that describes a polygon. 240 The last point is assumed ot join the first point. 241 This is in UTM (lat long would b better though) 242 243 ll_lat - lower left latitude, in decimal degrees 244 ll-long - lower left longitude, in decimal degrees 245 grid_spacing - in decimal degrees 246 247 """ 248 249 msg = "grid_spacing can not be zero" 250 assert not grid_spacing == 0, msg 251 252 a = boundary_polygon 253 254 # List of segments. Each segment is two points. 255 segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))] 256 257 # convert the segs to Lat's and longs. 258 # Don't assume the zone of the segments is the same as the lower left 259 # corner of the lat long data!! They can easily be in different zones 260 lat_long_set = frozenset() 261 for seg in segs: 262 points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing, 263 lat_amount, long_amount, zone, 264 isSouthHemisphere) 265 lat_long_set |= frozenset(points_lat_long) 266 267 if lat_long_set == frozenset([]): 268 msg = "URS region specified and polygon does not overlap." 269 raise ValueError, msg 270 271 # Warning there is no info in geospatial saying the hemisphere of 272 # these points. There should be. 273 geo = Geospatial_data(data_points=list(lat_long_set), 274 points_are_lats_longs=True) 275 276 return geo 277 278 279 ## 280 # @brief Get the points that are needed to interpolate any point a a segment. 281 # @param seg Two points in the UTM. 282 # @param ll_lat Lower left latitude, in decimal degrees 283 # @param ll_long Lower left longitude, in decimal degrees 284 # @param grid_spacing 285 # @param lat_amount 286 # @param long_amount 287 # @param zone 288 # @param isSouthHemisphere 289 # @return A list of points. 290 def points_needed(seg, ll_lat, ll_long, grid_spacing, 291 lat_amount, long_amount, zone, 292 isSouthHemisphere): 293 """ 294 seg is two points, in UTM 295 return a list of the points, in lats and longs that are needed to 296 interpolate any point on the segment. 297 """ 298 299 from math import sqrt 300 301 geo_reference = Geo_reference(zone=zone) 302 geo = Geospatial_data(seg, geo_reference=geo_reference) 303 seg_lat_long = geo.get_data_points(as_lat_long=True, 304 isSouthHemisphere=isSouthHemisphere) 305 306 # 1.415 = 2^0.5, rounded up.... 307 sqrt_2_rounded_up = 1.415 308 buffer = sqrt_2_rounded_up * grid_spacing 309 310 max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer 311 max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer 312 min_lat = min(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer 313 min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer 314 315 first_row = (min_long - ll_long) / grid_spacing 316 317 # To round up 318 first_row_long = int(round(first_row + 0.5)) 319 320 last_row = (max_long - ll_long) / grid_spacing # round down 321 last_row_long = int(round(last_row)) 322 323 first_row = (min_lat - ll_lat) / grid_spacing 324 # To round up 325 first_row_lat = int(round(first_row + 0.5)) 326 327 last_row = (max_lat - ll_lat) / grid_spacing # round down 328 last_row_lat = int(round(last_row)) 329 330 max_distance = 157147.4112 * grid_spacing 331 points_lat_long = [] 332 333 # Create a list of the lat long points to include. 334 for index_lat in range(first_row_lat, last_row_lat + 1): 335 for index_long in range(first_row_long, last_row_long + 1): 336 lat = ll_lat + index_lat*grid_spacing 337 long = ll_long + index_long*grid_spacing 338 339 #filter here to keep good points 340 if keep_point(lat, long, seg, max_distance): 341 points_lat_long.append((lat, long)) #must be hashable 342 343 # Now that we have these points, lets throw ones out that are too far away 344 return points_lat_long 345 346 347 348 ## 349 # @brief 350 # @param lat 351 # @param long 352 # @param seg Two points in UTM. 353 # @param max_distance 354 def keep_point(lat, long, seg, max_distance): 355 """ 356 seg is two points, UTM 357 """ 358 359 from math import sqrt 360 361 _ , x0, y0 = redfearn(lat, long) 362 x1 = seg[0][0] 363 y1 = seg[0][1] 364 x2 = seg[1][0] 365 y2 = seg[1][1] 366 x2_1 = x2-x1 367 y2_1 = y2-y1 368 try: 369 d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt( \ 370 (x2_1)*(x2_1)+(y2_1)*(y2_1)) 371 except ZeroDivisionError: 372 if sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1)) == 0 \ 373 and abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1)) == 0: 374 return True 375 else: 376 return False 377 378 return d <= max_distance 379 380 381
Note: See TracChangeset
for help on using the changeset viewer.