Changeset 7766

Timestamp:

Jun 2, 2010, 1:15:38 PM (15 years ago)

Author:

hudson

Message:

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

Location:

trunk/anuga_core/source/anuga/file

Files:

• mux.py (modified) (4 diffs)
• sts.py (modified) (1 diff)
• test_mux.py (modified) (3 diffs)
• test_urs.py (added)
• urs.py (modified) (2 diffs)

trunk/anuga_core/source/anuga/file/mux.py

@@ +1 @@
+"""
+    Read a mux2 file.
+"""
+
 from anuga.utilities.numerical_tools import ensure_numeric
 import numpy as num     
         
- ################################################################################
+################################################################################
 # READ MUX2 FILES line of points
 ################################################################################
+
+WAVEHEIGHT_MUX_LABEL = '-z-mux'
+EAST_VELOCITY_LABEL =  '-e-mux'
+NORTH_VELOCITY_LABEL =  '-n-mux'
 
 WAVEHEIGHT_MUX2_LABEL = '-z-mux2'
@@ -43 +51 @@
     # Convert verbose to int C flag
     if verbose:
-        verbose=1
+        verbose = 1
     else:
-        verbose=0
+        verbose = 0
 
     if weights is None:
@@ -102 +110 @@
     longitudes = num.zeros(number_of_selected_stations, num.float)
     elevation = num.zeros(number_of_selected_stations, num.float)
-    quantity = num.zeros((number_of_selected_stations, parameters_index), num.float)
+    quantity = num.zeros((number_of_selected_stations, parameters_index), \
+                                                    num.float)
 
     starttime = 1e16
@@ -116 +125 @@
 
 
-##
-# @brief ??
-# @param mux_times ??
-# @param mint ??
-# @param maxt ??
-# @return ??
-def read_time_from_mux(mux_times, mint, maxt):
-    """
-    """
 
-    if mint == None:
-        mux_times_start_i = 0
-    else:
-        mux_times_start_i = num.searchsorted(mux_times, mint)
-
-    if maxt == None:
-        mux_times_fin_i = len(mux_times)
-    else:
-        maxt += 0.5 # so if you specify a time where there is
-                    # data that time will be included
-        mux_times_fin_i = num.searchsorted(mux_times, maxt)
-
-    return mux_times_start_i, mux_times_fin_i
-
-
-       

trunk/anuga_core/source/anuga/file/sts.py

@@ -284 +284 @@
     except:
         msg = 'Cannot open %s' % stsname + '.sts'
-        raise msg
+        raise IOError(msg)
 
     xllcorner = fid.xllcorner[0]

trunk/anuga_core/source/anuga/file/test_mux.py

@@ -10 +10 @@
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 
+from mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
+                            NORTH_VELOCITY_LABEL
+
 from mux import WAVEHEIGHT_MUX2_LABEL, EAST_VELOCITY_MUX2_LABEL, \
                 NORTH_VELOCITY_MUX2_LABEL
@@ -15 +18 @@
 from mux import read_mux2_py
 from anuga.file_conversion.urs2sts import urs2sts
-
-class TestCase(unittest.TestCase):
+from anuga.file.urs import Read_urs
+
+class Test_Mux(unittest.TestCase):
     def setUp(self):
         pass
@@ -1487 +1491 @@
         
         self.delete_mux(files)
-        
+ 
+    def test_Urs_points(self):
+        time_step_count = 3
+        time_step = 2
+        lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,115)]
+        base_name, files = self.write_mux(lat_long_points,
+                                          time_step_count, time_step)
+        for file in files:
+            urs = Read_urs(file)
+            assert time_step_count == urs.time_step_count
+            assert time_step == urs.time_step
+
+            for lat_lon, dep in map(None, lat_long_points, urs.lonlatdep):
+                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
+                    assert num.allclose(n, dep[2])
+                        
+            count = 0
+            for slice in urs:
+                count += 1
+                #print slice
+                for lat_lon, quantity in map(None, lat_long_points, slice):
+                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
+                    #print "quantity", quantity
+                    #print "e", e
+                    #print "n", n
+                    if file[-5:] == WAVEHEIGHT_MUX_LABEL[-5:] or \
+                           file[-5:] == NORTH_VELOCITY_LABEL[-5:] :
+                        assert num.allclose(e, quantity)
+                    if file[-5:] == EAST_VELOCITY_LABEL[-5:]:
+                        assert num.allclose(n, quantity)
+            assert count == time_step_count
+                     
+        self.delete_mux(files)        
+
+
+
         
 ################################################################################
 
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Mux,'test')
     runner = unittest.TextTestRunner() #verbosity=2)
     runner.run(suite)

trunk/anuga_core/source/anuga/file/urs.py

@@ -1 @@
-##
-# @brief 
+from struct import pack, unpack
+import array as p_array
+import numpy as num
+
+
+from anuga.coordinate_transforms.geo_reference import Geo_reference
+
+from anuga.geospatial_data.geospatial_data import ensure_absolute, \
+                                                    Geospatial_data
+
+from anuga.coordinate_transforms.redfearn import redfearn
+
 class Read_urs:
     """
@@ -88 +98 @@
         self.mux_file.close()
    
+
+
+
+### PRODUCING THE POINTS NEEDED FILE ###
+
+# Ones used for FESA 2007 results
+#LL_LAT = -50.0
+#LL_LONG = 80.0
+#GRID_SPACING = 1.0/60.0
+#LAT_AMOUNT = 4800
+#LONG_AMOUNT = 3600
+
+
+##
+# @brief 
+# @param file_name 
+# @param boundary_polygon 
+# @param zone 
+# @param ll_lat 
+# @param ll_long 
+# @param grid_spacing 
+# @param lat_amount 
+# @param long_amount 
+# @param isSouthernHemisphere 
+# @param export_csv 
+# @param use_cache 
+# @param verbose True if this function is to be verbose.
+# @return 
+
+def save_boundary_as_urs(file_name, boundary_polygon, zone,
+                              ll_lat, ll_long,
+                              grid_spacing,
+                              lat_amount, long_amount,
+                              isSouthernHemisphere=True,
+                              export_csv=False, use_cache=False,
+                              verbose=False):
+    """
+    Given the info to replicate the URS grid and a polygon output
+    a file that specifies the cloud of boundary points for URS.
+
+    This creates a .urs file.  This is in the format used by URS;
+    1st line is the number of points,
+    each line after represents a point,in lats and longs.
+
+    Note: The polygon cannot cross zones or hemispheres.
+
+    A work-a-round for different zones or hemispheres is to run this twice,
+    once for each zone, and then combine the output.
+
+    file_name - name of the urs file produced for David.
+    boundary_polygon - a list of points that describes a polygon.
+                      The last point is assumed ot join the first point.
+                      This is in UTM (lat long would be better though)
+
+     This is info about the URS model that needs to be inputted.
+
+    ll_lat - lower left latitude, in decimal degrees
+    ll-long - lower left longitude, in decimal degrees
+    grid_spacing - in deciamal degrees
+    lat_amount - number of latitudes
+    long_amount- number of longs
+
+    Don't add the file extension.  It will be added.
+    """
+
+    geo = calculate_boundary_points(boundary_polygon, zone, ll_lat, ll_long,
+                            grid_spacing,
+                            lat_amount, long_amount, isSouthernHemisphere,
+                            use_cache, verbose)
+
+    if not file_name[-4:] == ".urs":
+        file_name += ".urs"
+
+    geo.export_points_file(file_name, isSouthHemisphere=isSouthernHemisphere)
+
+    if export_csv:
+        if file_name[-4:] == ".urs":
+            file_name = file_name[:-4] + ".csv"
+        geo.export_points_file(file_name)
+
+    return geo
+
+
+##
+# @brief 
+# @param boundary_polygon
+# @param zone
+# @param ll_lat
+# @param ll_long
+# @param grid_spacing
+# @param lat_amount
+# @param long_amount
+# @param isSouthHemisphere
+# @param use_cache
+# @param verbose
+def calculate_boundary_points(boundary_polygon, zone, ll_lat,
+                      ll_long, grid_spacing,
+                      lat_amount, long_amount, isSouthHemisphere=True,
+                      use_cache=False, verbose=False):
+    args = (boundary_polygon,
+            zone, ll_lat,
+            ll_long, grid_spacing,
+            lat_amount, long_amount, isSouthHemisphere)
+    kwargs = {}
+
+    if use_cache is True:
+        try:
+            from anuga.caching import cache
+        except:
+            msg = 'Caching was requested, but caching module' \
+                  'could not be imported'
+            raise msg
+
+        geo = cache(_calculate_boundary_points,
+                    args, kwargs,
+                    verbose=verbose,
+                    compression=False)
+    else:
+        geo = apply(_calculate_boundary_points, args, kwargs)
+
+    return geo
+
+
+##
+# @brief 
+# @param boundary_polygon An iterable of points that describe a polygon.
+# @param zone
+# @param ll_lat Lower left latitude, in decimal degrees
+# @param ll_long Lower left longitude, in decimal degrees
+# @param grid_spacing Grid spacing in decimal degrees.
+# @param lat_amount
+# @param long_amount
+# @param isSouthHemisphere
+def _calculate_boundary_points(boundary_polygon,
+                       zone, ll_lat,
+                       ll_long, grid_spacing,
+                       lat_amount, long_amount,
+                       isSouthHemisphere):
+    """
+    boundary_polygon - a list of points that describes a polygon.
+                      The last point is assumed ot join the first point.
+                      This is in UTM (lat long would b better though)
+
+    ll_lat - lower left latitude, in decimal degrees
+    ll-long - lower left longitude, in decimal degrees
+    grid_spacing - in decimal degrees
+
+    """
+
+    msg = "grid_spacing can not be zero"
+    assert not grid_spacing == 0, msg
+
+    a = boundary_polygon
+
+    # List of segments.  Each segment is two points.
+    segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))]
+
+    # convert the segs to Lat's and longs.
+    # Don't assume the zone of the segments is the same as the lower left
+    # corner of the lat long data!!  They can easily be in different zones
+    lat_long_set = frozenset()
+    for seg in segs:
+        points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing,
+                                        lat_amount, long_amount, zone,
+                                        isSouthHemisphere)
+        lat_long_set |= frozenset(points_lat_long)
+
+    if lat_long_set == frozenset([]):
+        msg = "URS region specified and polygon does not overlap."
+        raise ValueError, msg
+
+    # Warning there is no info in geospatial saying the hemisphere of
+    # these points.  There should be.
+    geo = Geospatial_data(data_points=list(lat_long_set),
+                          points_are_lats_longs=True)
+
+    return geo
+
+
+##
+# @brief Get the points that are needed to interpolate any point a a segment.
+# @param seg Two points in the UTM.
+# @param ll_lat Lower left latitude, in decimal degrees
+# @param ll_long Lower left longitude, in decimal degrees
+# @param grid_spacing 
+# @param lat_amount 
+# @param long_amount 
+# @param zone 
+# @param isSouthHemisphere 
+# @return A list of points.
+def points_needed(seg, ll_lat, ll_long, grid_spacing,
+                  lat_amount, long_amount, zone,
+                  isSouthHemisphere):
+    """
+    seg is two points, in UTM
+    return a list of the points, in lats and longs that are needed to
+    interpolate any point on the segment.
+    """
+
+    from math import sqrt
+
+    geo_reference = Geo_reference(zone=zone)
+    geo = Geospatial_data(seg, geo_reference=geo_reference)
+    seg_lat_long = geo.get_data_points(as_lat_long=True,
+                                       isSouthHemisphere=isSouthHemisphere)
+
+    # 1.415 = 2^0.5, rounded up....
+    sqrt_2_rounded_up = 1.415
+    buffer = sqrt_2_rounded_up * grid_spacing
+
+    max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer
+    max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer
+    min_lat = min(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer
+    min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
+
+    first_row = (min_long - ll_long) / grid_spacing
+
+    # To round up
+    first_row_long = int(round(first_row + 0.5))
+
+    last_row = (max_long - ll_long) / grid_spacing # round down
+    last_row_long = int(round(last_row))
+
+    first_row = (min_lat - ll_lat) / grid_spacing
+    # To round up
+    first_row_lat = int(round(first_row + 0.5))
+
+    last_row = (max_lat - ll_lat) / grid_spacing # round down
+    last_row_lat = int(round(last_row))
+
+    max_distance = 157147.4112 * grid_spacing
+    points_lat_long = []
+
+    # Create a list of the lat long points to include.
+    for index_lat in range(first_row_lat, last_row_lat + 1):
+        for index_long in range(first_row_long, last_row_long + 1):
+            lat = ll_lat + index_lat*grid_spacing
+            long = ll_long + index_long*grid_spacing
+
+            #filter here to keep good points
+            if keep_point(lat, long, seg, max_distance):
+                points_lat_long.append((lat, long)) #must be hashable
+
+    # Now that we have these points, lets throw ones out that are too far away
+    return points_lat_long
+       
+
+
+##
+# @brief 
+# @param lat
+# @param long
+# @param seg Two points in UTM.
+# @param max_distance
+def keep_point(lat, long, seg, max_distance):
+    """
+    seg is two points, UTM
+    """
+
+    from math import sqrt
+
+    _ , x0, y0 = redfearn(lat, long)
+    x1 = seg[0][0]
+    y1 = seg[0][1]
+    x2 = seg[1][0]
+    y2 = seg[1][1]
+    x2_1 = x2-x1
+    y2_1 = y2-y1
+    try:
+        d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt( \
+            (x2_1)*(x2_1)+(y2_1)*(y2_1))
+    except ZeroDivisionError:
+        if sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1)) == 0 \
+           and abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1)) == 0:
+            return True
+        else:
+            return False
+
+    return d <= max_distance
+
+
+