Changeset 4223


Ignore:
Timestamp:
Feb 6, 2007, 2:28:47 PM (17 years ago)
Author:
duncan
Message:

Start of improved urs2sww, along with other stuff.

Location:
anuga_core/source/anuga
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/coordinate_transforms/redfearn.py

    r4202 r4223  
    189189
    190190    Assume the false_easting and false_northing are the same for each list.
    191     If points end up belonging to different UTM zones, an ANUGAerror is thrown.   
     191    If points end up in different UTM zones, an ANUGAerror is thrown.   
    192192    """
    193193
  • anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py

    r4211 r4223  
    6868              use_least_squares=False,
    6969              use_file_type=None,
    70               save=False):
     70              save=False,
     71              verbose=False):
    7172        '''
    7273        num_of_points
     
    117118                                         points_dict['point_attributes'])
    118119                    G1.export_points_file(fileName, absolute=True)
    119                     calc = interp.fit(fileName)
     120                    calc = interp.fit(fileName), verbose=verbose)
    120121                    os.remove(fileName)
    121122                   
  • anuga_core/source/anuga/fit_interpolate/run_long_benchmark.py

    r4211 r4223  
    1919#num_of_points_list = [10]
    2020#maxArea_list = [0.1, 0.001]
    21 num_of_points_list = [10, 100] #, 10000, 100000] #, 10000000]
    22 maxArea_list = [0.1, 0.001] #, 0.00001, 0.0000001]
     21num_of_points_list = [10000, 100000] #, 10000, 100000] #, 10000000]
     22maxArea_list = [0.1, 0.01, 0.001] #, 0.00001, 0.0000001]
    2323max_points_per_cell_list = [8]
    24 use_file_type_list = [None,'txt']
     24use_file_type_list = [None,'txt', 'pts']
    2525
    2626fd = open(ofile,'a')
     
    3737
    3838
    39 for use_file_type in use_file_type_list:
    40     for is_fit in is_fit_list:
    41         for num_of_points in num_of_points_list:
    42             for maxArea in maxArea_list:
     39for maxArea in maxArea_list:
     40    for use_file_type in use_file_type_list:
     41        for is_fit in is_fit_list:
     42            for num_of_points in num_of_points_list:
    4343                for max_points_per_cell in max_points_per_cell_list:
    4444   
     
    5252                    print "mem", mem
    5353                    fd.write(str(use_file_type) + delimiter +
    54                              str(is_fit) + delimiter +
    5554                             str(num_of_points) + delimiter +
    5655                             str(maxArea) + delimiter +
    5756                             str(num_tri) + delimiter +
    5857                             str(max_points_per_cell) + delimiter +
     58                             str(is_fit) + delimiter +
    5959                             str(mem)  + delimiter +
    6060                             str(time) + delimiter + "\n")
  • anuga_core/source/anuga/geospatial_data/geospatial_data.py

    r4202 r4223  
    1414
    1515from Scientific.IO.NetCDF import NetCDFFile   
    16    
     16from anuga.coordinate_transforms.lat_long_UTM_conversion import UTMtoLL   
    1717from anuga.utilities.numerical_tools import ensure_numeric
    18 from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
     18from anuga.coordinate_transforms.geo_reference import Geo_reference, \
     19     TitleError, DEFAULT_ZONE
    1920from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
    2021from anuga.utilities.anuga_exceptions import ANUGAError
     
    332333        return self.geo_reference
    333334       
    334     def get_data_points(self, absolute=True, geo_reference=None):
     335    def get_data_points(self, absolute=True, geo_reference=None,
     336                        as_lat_long=False):
    335337        """Get coordinates for all data points as an Nx2 array
    336338
     
    344346        Default: absolute is True.
    345347        """
    346 
     348        if as_lat_long is True:
     349            msg = "Points need a zone to be converted into lats and longs"
     350            assert self.geo_reference is not None, msg
     351            zone = self.geo_reference.get_zone()
     352            assert self.geo_reference.get_zone() is not DEFAULT_ZONE, msg
     353            lats_longs = []
     354            for point in self.get_data_points(True):
     355                ### UTMtoLL(northing, easting, zone,
     356                lat_calced, long_calced = UTMtoLL(point[1],point[0], zone)
     357                lats_longs.append([lat_calced, long_calced])
     358            return lats_longs
     359           
    347360        if absolute is True and geo_reference is None:
    348361            return self.geo_reference.get_absolute(self.data_points)
  • anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

    r4202 r4223  
    77from math import sqrt, pi
    88import tempfile
     9from sets import ImmutableSet
    910
    1011from anuga.geospatial_data.geospatial_data import *
     
    162163        assert allclose(results, points_rel)
    163164
    164        
     165 
     166    def test_get_data_points_lat_long(self):
     167        # lat long [-30.],[130]
     168        #Zone:   52   
     169        #Easting:  596450.153  Northing: 6680793.777
     170        # lat long [-32.],[131]
     171        #Zone:   52   
     172        #Easting:  688927.638  Northing: 6457816.509
     173       
     174        points_Lat_long = [[-30.,130], [-32,131]]
     175       
     176        spatial = Geospatial_data(latitudes=[-30, -32.],
     177                                  longitudes=[130, 131])
     178
     179        results = spatial.get_data_points(as_lat_long=True)
     180        #print "test_get_data_points_lat_long - results", results
     181        #print "points_Lat_long",points_Lat_long
     182        assert allclose(results, points_Lat_long)
     183     
     184    def test_get_data_points_lat_longII(self):
     185        # x,y  North,east long,lat
     186        boundary_polygon = [[ 250000, 7630000]]
     187        zone = 50
     188       
     189        geo_reference = Geo_reference(zone=zone)
     190        geo = Geospatial_data(boundary_polygon,geo_reference=geo_reference)
     191        seg_lat_long = geo.get_data_points(as_lat_long=True)
     192        lat_result = degminsec2decimal_degrees(-21,24,54)
     193        long_result = degminsec2decimal_degrees(114,35,17.89)
     194        #print "seg_lat_long", seg_lat_long [0][0]
     195        #print "lat_result",lat_result
     196        assert allclose(seg_lat_long[0][0], lat_result)#lat
     197        assert allclose(seg_lat_long[0][1], long_result)#long
     198
     199             
    165200    def test_set_geo_reference(self):
    166201        points_ab = [[12.5,34.7],[-4.5,-60.0]]
     
    10581093
    10591094    def verbose_test_load_pts_blocking(self):
    1060         """ test_load_csv(self):
    1061         space delimited
    1062         """
     1095       
    10631096        import os
    10641097       
     
    11301163                        [10.0, 0.0])
    11311164        assert allclose(geo_list[1].get_data_points(),
    1132                         [[1.0, 0.0]])       
     1165                        [[1.0, 0.0],[0.0, 1.0] ])       
    11331166        assert allclose(geo_list[1].get_attributes(attribute_name='elevation'),
    1134                         [10.4])
     1167                        [10.0, 0.0])
    11351168           
    11361169        os.remove(fileName) 
     
    19812014            self.failUnless(0 ==1,  'Error not thrown error!')
    19822015
     2016
     2017    def test_lat_long_set(self):
     2018        lat_gong = degminsec2decimal_degrees(-34,30,0.)
     2019        lon_gong = degminsec2decimal_degrees(150,55,0.)
     2020       
     2021        lat_2 = degminsec2decimal_degrees(-34,00,0.)
     2022        lon_2 = degminsec2decimal_degrees(150,00,0.)
     2023        p1 = (lat_gong, lon_gong)
     2024        p2 = (lat_2, lon_2)
     2025        points = ImmutableSet([p1, p2, p1])
     2026        gsd = Geospatial_data(data_points=list(points),
     2027                              points_are_lats_longs=True)
     2028
     2029        points = gsd.get_data_points(absolute=True)
     2030       
     2031        assert allclose(points[0][0], 308728.009)
     2032        assert allclose(points[0][1], 6180432.601)
     2033        assert allclose(points[1][0],  222908.705)
     2034        assert allclose(points[1][1], 6233785.284)
     2035        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
     2036                        'Bad zone error!')
     2037        points = gsd.get_data_points(as_lat_long=True)
     2038        #print "test_lat_long_set points", points
     2039        assert allclose(points[1][0], -34)
     2040        assert allclose(points[1][1], 150)
     2041
    19832042    def test_len(self):
    19842043       
     
    20252084if __name__ == "__main__":
    20262085
    2027 #    suite = unittest.makeSuite(Test_Geospatial_data, 'test_split')
    2028 #    suite = unittest.makeSuite(Test_Geospatial_data, 'test_clip0')
    2029     suite = unittest.makeSuite(Test_Geospatial_data, 'test')
     2086    suite = unittest.makeSuite(Test_Geospatial_data, 'test_lat_long_set')
     2087    #suite = unittest.makeSuite(Test_Geospatial_data, 'verbose_test_load_pts_blocking')
     2088    #suite = unittest.makeSuite(Test_Geospatial_data, 'test')
    20302089    runner = unittest.TextTestRunner()
    20312090    runner.run(suite)
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r4059 r4223  
    6868
    6969
     70from anuga.coordinate_transforms.redfearn import redfearn
    7071from anuga.coordinate_transforms.geo_reference import Geo_reference
    7172from anuga.geospatial_data.geospatial_data import Geospatial_data
     
    44044405    return long, lat, quantity
    44054406
    4406     ####  END URS 2 SWW  ###     
     4407    ####  END URS 2 SWW  ###
     4408
     4409    #### URS UNGRIDDED 2 SWW ###
     4410
     4411def URS_points_needed_to_file(boundary_polygon, ll_lat, ll_long,
     4412                              grid_spacing,
     4413                      lat_amount, long_amount, zone=None):
     4414   
     4415    geo = URS_points_needed(boundary_polygon, ll_lat, ll_long, grid_spacing,
     4416                      lat_amount, long_amount, zone=zone)
     4417   
     4418def URS_points_needed(boundary_polygon, ll_lat, ll_long, grid_spacing,
     4419                      lat_amount, long_amount, zone=None):
     4420    """
     4421
     4422    boundary_polygon - a list of points that describes a polygon.
     4423                      The last point is assumed ot join the first point.
     4424                      This is in UTM (lat long would be better though)
     4425
     4426    ll_lat - lower left latitude, in decimal degrees
     4427    ll-long - lower left longitude, in decimal degrees
     4428    grid_spacing - in deciamal degrees
     4429
     4430    """
     4431    from sets import ImmutableSet
     4432   
     4433    msg = "grid_spacing can not be zero"
     4434    assert not grid_spacing ==0, msg
     4435    a = boundary_polygon
     4436    # List of segments.  Each segment is two points.
     4437    segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))]
     4438
     4439    # convert the segs to Lat's and longs.
     4440    if zone is None:
     4441        # Assume the zone of the segments is the same as the lower left
     4442        # corner of the lat long data
     4443        zone, _, _ = redfearn(ll_lat, ll_long)
     4444    lat_long_set = ImmutableSet()
     4445    for seg in segs:
     4446        points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing,
     4447                      lat_amount, long_amount, zone)
     4448        lat_long_set |= ImmutableSet(points_lat_long)
     4449    #print "lat_long_set",lat_long_set
     4450    geo = Geospatial_data(data_points=list(lat_long_set),
     4451                              points_are_lats_longs=True)
     4452    return geo
     4453   
     4454def points_needed(seg, ll_lat, ll_long, grid_spacing,
     4455                  lat_amount, long_amount, zone):
     4456    """
     4457    return a list of the points, in lats and longs that are needed to
     4458    interpolate any point on the segment.
     4459    """
     4460    #print "zone",zone
     4461    geo_reference = Geo_reference(zone=zone)
     4462    #print "seg",seg
     4463    geo = Geospatial_data(seg,geo_reference=geo_reference)
     4464    seg_lat_long = geo.get_data_points(as_lat_long=True)
     4465    #print "seg_lat_long", seg_lat_long
     4466    # 1.415 = 2^0.5, rounded up....
     4467    buffer = 1.415 * grid_spacing
     4468   
     4469    #
     4470   
     4471    max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer
     4472    max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer
     4473    min_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer
     4474    min_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
     4475
     4476    #print "ll_lat", ll_lat
     4477    #print "ll_long", ll_long
     4478    #print "max_lat", max_lat
     4479    #print "max_long", max_long
     4480    #print "min_lat", min_lat
     4481    #print "min_long", min_long
     4482    first_row = (min_long - ll_long)/grid_spacing
     4483    # To round up
     4484    first_row_long = int(round(first_row + 0.5))
     4485    #print "first_row", first_row_long
     4486
     4487    last_row = (max_long - ll_long)/grid_spacing # round down
     4488    last_row_long = int(round(last_row))
     4489    #print "last_row",last_row _long
     4490
     4491   
     4492    first_row = (min_lat - ll_lat)/grid_spacing
     4493    # To round up
     4494    first_row_lat = int(round(first_row + 0.5))
     4495    #print "first_row", first_row_lat
     4496
     4497    last_row = (max_lat - ll_lat)/grid_spacing # round down
     4498    last_row_lat = int(round(last_row))
     4499    #print "last_row",last_row_lat
     4500
     4501    points_lat_long = []
     4502    # Create a list of the lat long points to include.
     4503    for index_lat in range(first_row_lat, last_row_lat + 1):
     4504        for index_long in range(first_row_long, last_row_long + 1):
     4505            lat = ll_lat + index_lat*grid_spacing
     4506            long = ll_long + index_long*grid_spacing
     4507            points_lat_long.append((lat, long)) #must be hashable
     4508    #print "points_lat_long", points_lat_long
     4509    return points_lat_long   
     4510   
     4511    #### END URS UNGRIDDED 2 SWW ###
     4512
    44074513#-------------------------------------------------------------
    44084514if __name__ == "__main__":
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r4202 r4223  
    50355035            raise msg
    50365036       
    5037     def trial_loading(self):
    5038         basename_in = 'karratha'
    5039         basename_out = basename_in
    5040         urs2sww(basename_in, basename_out, remove_nc_files=True,
    5041                 zscale=10000000)
    50425037    #### END TESTS FOR URS 2 SWW  ###
    50435038
     5039    #### TESTS URS UNGRIDDED 2 SWW ###
     5040    def test_URS_points_needed(self):
     5041        ll_lat = -21.5
     5042        ll_long = 114.5
     5043        grid_spacing = 1./60.
     5044        lat_amount = 30
     5045        long_amount = 30
     5046
     5047        boundary_polygon = [[250000,7660000],[280000,7660000],
     5048                             [280000,7630000],[250000,7630000]]
     5049        geo=URS_points_needed(boundary_polygon,
     5050                              ll_lat, ll_long, grid_spacing,
     5051                              lat_amount, long_amount)
     5052        # to test this geo, can info from it be transfered onto the boundary
     5053        # poly?
     5054        #Maybe see if I can fit the data to the polygon - have to represent
     5055        # the poly as points though.
     5056       
     5057        results = geo.get_data_points(as_lat_long=True)
     5058        #print 'results',results
     5059       
     5060    def X_test_URS_points_neededII(self):
     5061        ll_lat = -21.5
     5062        ll_long = 114.5
     5063        grid_spacing = 1./60.
     5064        lat_amount = 30
     5065        long_amount = 30
     5066
     5067        # change this so lats and longs are inputed, then converted
     5068        #boundary_polygon = [[7660000,250000],[7660000,280000],
     5069        #                     [7630000,280000],[7630000,250000]]
     5070        URS_points_needed(boundary_polygon, ll_lat, ll_long, grid_spacing, \
     5071                      lat_amount, long_amount)
     5072       
     5073    #### END TESTS URS UNGRIDDED 2 SWW ###
    50445074       
    50455075#-------------------------------------------------------------
    50465076if __name__ == "__main__":
    5047     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sww_m')
     5077    #suite = unittest.makeSuite(Test_Data_Manager,'test_URS_points_needed')
    50485078    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_min_max_indexes_lat_ascending')
    50495079    #suite = unittest.makeSuite(Test_Data_Manager,'test_ferret2sww_lat_long')
Note: See TracChangeset for help on using the changeset viewer.