source: trunk/anuga_core/source/anuga/file/urs.py @ 7859

Last change on this file since 7859 was 7766, checked in by hudson, 14 years ago

File tests passing OK - extra module for .urs files.

File size: 11.5 KB
Line 
1from struct import pack, unpack
2import array as p_array
3import numpy as num
4
5
6from anuga.coordinate_transforms.geo_reference import Geo_reference
7
8from anuga.geospatial_data.geospatial_data import ensure_absolute, \
9                                                    Geospatial_data
10
11from anuga.coordinate_transforms.redfearn import redfearn
12
13class Read_urs:
14    """
15    Read the info in URS mux files.
16
17    for the quantities here's a correlation between the file names and
18    what they mean;
19    z-mux is height above sea level, m
20    e-mux is velocity is Eastern direction, m/s
21    n-mux is velocity is Northern direction, m/s
22    """
23
24    ##
25    # @brief Initialize this instance of Urs_points.
26    # @param urs_file Path to the underlying data file.
27    def __init__(self, urs_file):
28        self.iterated = False
29        columns = 3                         # long, lat , depth
30        mux_file = open(urs_file, 'rb')
31
32        # Number of points/stations
33        (self.points_num,) = unpack('i', mux_file.read(4))
34
35        # nt, int - Number of time steps
36        (self.time_step_count,) = unpack('i', mux_file.read(4))
37        #dt, float - time step, seconds
38        (self.time_step,) = unpack('f', mux_file.read(4))
39        msg = "Bad data in the urs file."
40        if self.points_num < 0:
41            mux_file.close()
42            raise ANUGAError, msg
43        if self.time_step_count < 0:
44            mux_file.close()
45            raise ANUGAError, msg
46        if self.time_step < 0:
47            mux_file.close()
48            raise ANUGAError, msg
49
50        # The depth is in meters, and it is the distance from the ocean
51        # to the sea bottom.
52        lonlatdep = p_array.array('f')
53        lonlatdep.read(mux_file, columns * self.points_num)
54        lonlatdep = num.array(lonlatdep, dtype=num.float)
55        lonlatdep = num.reshape(lonlatdep, (self.points_num, columns))
56        self.lonlatdep = lonlatdep
57
58        self.mux_file = mux_file
59        # check this array
60
61    ##
62    # @brief Allow iteration over quantity data wrt time.
63    def __iter__(self):
64        """
65        iterate over quantity data which is with respect to time.
66
67        Note: You can only iterate once over an object
68
69        returns quantity infomation for each time slice
70        """
71
72        msg =  "You can only interate once over a urs file."
73        assert not self.iterated, msg
74
75        self.iter_time_step = 0
76        self.iterated = True
77
78        return self
79
80    ##
81    # @brief
82    def next(self):
83        if self.time_step_count == self.iter_time_step:
84            self.close()
85            raise StopIteration
86
87        #Read in a time slice from mux file
88        hz_p_array = p_array.array('f')
89        hz_p_array.read(self.mux_file, self.points_num)
90        hz_p = num.array(hz_p_array, dtype=num.float)
91        self.iter_time_step += 1
92
93        return hz_p
94
95    ##
96    # @brief Close the mux file.
97    def close(self):
98        self.mux_file.close()
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
129def 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
195def 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
233def _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.
290def 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
354def 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 TracBrowser for help on using the repository browser.