Changeset 7766


Ignore:
Timestamp:
Jun 2, 2010, 1:15:38 PM (14 years ago)
Author:
hudson
Message:

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

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
    15from anuga.utilities.numerical_tools import ensure_numeric
    26import numpy as num     
    37       
    4  ################################################################################
     8################################################################################
    59# READ MUX2 FILES line of points
    610################################################################################
     11
     12WAVEHEIGHT_MUX_LABEL = '-z-mux'
     13EAST_VELOCITY_LABEL =  '-e-mux'
     14NORTH_VELOCITY_LABEL =  '-n-mux'
    715
    816WAVEHEIGHT_MUX2_LABEL = '-z-mux2'
     
    4351    # Convert verbose to int C flag
    4452    if verbose:
    45         verbose=1
     53        verbose = 1
    4654    else:
    47         verbose=0
     55        verbose = 0
    4856
    4957    if weights is None:
     
    102110    longitudes = num.zeros(number_of_selected_stations, num.float)
    103111    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)
    105114
    106115    starttime = 1e16
     
    116125
    117126
    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     """
    127127
    128     if mint == None:
    129         mux_times_start_i = 0
    130     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 is
    137                     # data that time will be included
    138         mux_times_fin_i = num.searchsorted(mux_times, maxt)
    139 
    140     return mux_times_start_i, mux_times_fin_i
    141 
    142 
    143        
  • trunk/anuga_core/source/anuga/file/sts.py

    r7762 r7766  
    284284    except:
    285285        msg = 'Cannot open %s' % stsname + '.sts'
    286         raise msg
     286        raise IOError(msg)
    287287
    288288    xllcorner = fid.xllcorner[0]
  • trunk/anuga_core/source/anuga/file/test_mux.py

    r7762 r7766  
    1010from anuga.coordinate_transforms.geo_reference import Geo_reference
    1111
     12from mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
     13                            NORTH_VELOCITY_LABEL
     14
    1215from mux import WAVEHEIGHT_MUX2_LABEL, EAST_VELOCITY_MUX2_LABEL, \
    1316                NORTH_VELOCITY_MUX2_LABEL
     
    1518from mux import read_mux2_py
    1619from anuga.file_conversion.urs2sts import urs2sts
    17 
    18 class TestCase(unittest.TestCase):
     20from anuga.file.urs import Read_urs
     21
     22class Test_Mux(unittest.TestCase):
    1923    def setUp(self):
    2024        pass
     
    14871491       
    14881492        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
    14901529       
    14911530################################################################################
    14921531
    14931532if __name__ == "__main__":
    1494     suite = unittest.makeSuite(TestCase,'test')
     1533    suite = unittest.makeSuite(Test_Mux,'test')
    14951534    runner = unittest.TextTestRunner() #verbosity=2)
    14961535    runner.run(suite)
  • trunk/anuga_core/source/anuga/file/urs.py

    r7760 r7766  
    1 ##
    2 # @brief
     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
    313class Read_urs:
    414    """
     
    8898        self.mux_file.close()
    8999   
     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 TracChangeset for help on using the changeset viewer.