Changeset 5361


Ignore:
Timestamp:
May 24, 2008, 11:55:52 AM (17 years ago)
Author:
jakeman
Message:

Adding ability to read sts file and create file function

Location:
anuga_core/source/anuga
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r5349 r5361  
    1919from anuga.geospatial_data.geospatial_data import ensure_absolute
    2020from math import sqrt, atan, degrees
    21 from exceptions import IOError
     21
    2222
    2323
    … …  
    110110    # Instead we pass in those attributes that are needed (and return them
    111111    # if modified)
    112 
    113112    kwargs = {'quantities': quantities,
    114113              'interpolation_points': interpolation_points,
    … …  
    179178    try:
    180179        fid = open(filename)
    181     except IOError, e:
     180    except Exception, e:
    182181        msg = 'File "%s" could not be opened: Error="%s"'\
    183182                  %(filename, e)
    184         raise IOError, msg # So IOErrors can be caught
     183        raise msg
    185184
    186185    line = fid.readline()
    … …  
    244243        raise Exception, msg
    245244
    246 
     245 
    247246    if interpolation_points is not None:
    248247        interpolation_points = ensure_absolute(interpolation_points)
    … …  
    250249        assert interpolation_points.shape[1] == 2, msg
    251250
    252     if verbose:
    253         print 'File_function using quantities %s from file %s' %(str(quantity_names), filename)
    254251
    255252    # Now assert that requested quantitites (and the independent ones)
    … …  
    273270
    274271    if filename[-3:] == 'tms' and spatial is True:
    275         msg = 'Files of type tms must not contain spatial information'
     272        msg = 'Files of type tms must not contain spatial  information'
    276273        raise msg
    277274
    278275    if filename[-3:] == 'sww' and spatial is False:
    279276        msg = 'Files of type sww must contain spatial information'       
    280         raise msg       
     277        raise msg
     278
     279    if filename[-3:] == 'sts' and spatial is False:
     280        #What if mux file only contains one point
     281        msg = 'Files of type sts must contain spatial information'       
     282        raise msg
    281283
    282284    # Get first timestep
    … …  
    300302        x = fid.variables['x'][:]
    301303        y = fid.variables['y'][:]
    302         triangles = fid.variables['volumes'][:]
     304        if filename[-3:] == 'sww':
     305            triangles = fid.variables['volumes'][:]
    303306
    304307        x = reshape(x, (len(x),1))
    … …  
    308311        if interpolation_points is not None:
    309312            # Adjust for georef
    310 
    311             # FIXME (Ole): Use geo_reference.get_relative
    312313            interpolation_points[:,0] -= xllcorner
    313314            interpolation_points[:,1] -= yllcorner       
    314        
    315 
    316 
    317315
    318316    if domain_starttime is not None:
    … …  
    350348
    351349    if not spatial:
    352         vertex_coordinates = triangles = interpolation_points = None         
    353 
     350        vertex_coordinates = triangles = interpolation_points = None
     351    if filename[-3:] == 'sts':#added
     352        triangles = None
     353        #vertex coordinates is position of urs gauges
    354354
    355355    # Return Interpolation_function instance as well as
    … …  
    782782        #print 'label', label
    783783        leg_label.append(label)
    784 
    785784
    786785       
    … …  
    21452144       
    21462145        callable_sww = file_function(sww_file,
    2147                                      quantities=core_quantities,
    2148                                      interpolation_points=points_array,
    2149                                      verbose=verbose,
    2150                                      use_cache=use_cache)
     2146                                 quantities=core_quantities,
     2147                                 interpolation_points=points_array,
     2148                                 verbose=verbose,
     2149                                 use_cache=use_cache)
    21512150
    21522151    gauge_file='gauge_'

  • anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

    r5321 r5361  
    9292            max_vertices_per_cell = MAX_VERTICES_PER_CELL
    9393        if mesh is None:
    94             # Fixme (DSG) Throw errors if triangles or vertex_coordinates
    95             # are None
     94            if vertex_coordinates is not None and  triangles is not None:
     95                # Fixme (DSG) Throw errors if triangles or vertex_coordinates
     96                # are None
    9697           
    97             #Convert input to Numeric arrays
    98             triangles = ensure_numeric(triangles, Int)
    99             vertex_coordinates = ensure_absolute(vertex_coordinates,
     98                #Convert input to Numeric arrays
     99                triangles = ensure_numeric(triangles, Int)
     100                vertex_coordinates = ensure_absolute(vertex_coordinates,
    100101                                                 geo_reference = mesh_origin)
    101102
    102             if verbose: print 'FitInterpolate: Building mesh'       
    103             self.mesh = Mesh(vertex_coordinates, triangles)
    104             #self.mesh.check_integrity() # Time consuming
     103                if verbose: print 'FitInterpolate: Building mesh'       
     104                self.mesh = Mesh(vertex_coordinates, triangles)
     105                #self.mesh.check_integrity() # Time consuming
     106            else:
     107                self.mesh = None
    105108        else:
    106109            self.mesh = mesh
     110
     111        if self.mesh is not None:
     112            if verbose: print 'FitInterpolate: Building quad tree'
     113            #This stores indices of vertices
     114            t0 = time.time()
     115            #print "self.mesh.get_extent(absolute=True)", \
     116            #self.mesh.get_extent(absolute=True)
     117            self.root = build_quadtree(self.mesh,
     118                                       max_points_per_cell = max_vertices_per_cell)
     119            #print "self.root",self.root.show()
    107120       
    108         if verbose: print 'FitInterpolate: Building quad tree'
    109         # This stores indices of vertices
    110         t0 = time.time()
    111         #print "self.mesh.get_extent(absolute=True)", \
    112               #self.mesh.get_extent(absolute=True)
    113         self.root = build_quadtree(self.mesh,
    114                                    max_points_per_cell = max_vertices_per_cell)
    115         #print "self.root",self.root.show()
    116        
    117         build_quadtree_time =  time.time()-t0
    118         set_last_triangle()
     121            build_quadtree_time =  time.time()-t0
     122            set_last_triangle()
    119123       
    120124    def __repr__(self):

  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r5321 r5361  
    8989                                verbose=verbose,
    9090                                max_vertices_per_cell=max_vertices_per_cell)
     91
     92    def interpolate_polyline(self,
     93                             f,
     94                             vertex_coordinates,
     95                             point_coordinates=None,
     96                             start_blocking_len=500000,
     97                             verbose=False):
     98 
     99        if isinstance(point_coordinates, Geospatial_data):
     100            point_coordinates = point_coordinates.get_data_points( \
     101                absolute = True)
     102 
     103        from utilities.polygon import point_on_line,point_on_line_py
     104        from Numeric import ones
     105        z=ones(len(point_coordinates),Float)
     106
     107        msg='point coordinates are not given (interpolate.py)'
     108        assert point_coordinates is not None, msg
     109        msg='function value must be specified at every interpolation node'
     110        assert f.shape[0]==vertex_coordinates.shape[0],msg
     111        msg='Must define function value at one or more nodes'
     112        assert f.shape[0]>0,msg
     113
     114        #print point_on_line_py(point_coordinates[3],[vertex_coordinates[0],vertex_coordinates[1]])
     115        #print vertex_coordinates[0],vertex_coordinates[1]
     116
     117        #print point_coordinates
     118        #print vertex_coordinates
     119        n=f.shape[0]
     120        if n==1:
     121            z=f*z
     122        elif n>1:
     123            index=0#index of segment on which last point was found
     124            k=0#number of points on line
     125            for i in range(len(point_coordinates)):
     126                found = False
     127                #find the segment the cetnroid lies on
     128                #Start with the segment the last point was found on
     129                #For n points there will be n-1 segments
     130                for j in range(n-1):
     131                    #print 'searcing segment', index+j
     132                    #reached last segment look at first segment
     133                    if (index+j)==n-2:
     134                        index=0
     135                    if point_on_line_py(point_coordinates[i],[vertex_coordinates[(j)+index],vertex_coordinates[(j+1)+index]]):
     136                        found=True
     137                        x0=vertex_coordinates[j][0];y0=vertex_coordinates[j][1]
     138                        x1=vertex_coordinates[j+1][0];y1=vertex_coordinates[j+1][1]
     139                        x2=point_coordinates[i][0];y2=point_coordinates[i][1]
     140                       
     141                        segment_len=sqrt((x1-x0)**2+(y1-y0)**2)
     142                        dist=sqrt((x2-x0)**2+(y2-y0)**2)
     143                        z[i]=(f[j+1]-f[j])/segment_len*dist+f[j]
     144                        #print 'element found on segment',j+index
     145                        index=j
     146                        break
     147                   
     148               
     149                if not found:
     150                    z[i]=0.0
     151                    #print 'point not on urs boundary'
     152                else:
     153                    k+=1
     154        #print 'number of midpoints on urs boundary',k
     155        #print z
     156        return z
    91157
    92158    # FIXME: What is a good start_blocking_len value?
    … …  
    529595            vertex_coordinates = ensure_absolute(vertex_coordinates)
    530596
    531             assert triangles is not None, 'Triangles array must be specified'
    532             triangles = ensure_numeric(triangles)
    533             self.spatial = True           
     597            if triangles is not None:
     598                triangles = ensure_numeric(triangles)
     599            self.spatial = True         
    534600
    535601        # Thin timesteps if needed
    … …  
    557623        # Precomputed spatial interpolation if requested
    558624        if interpolation_points is not None:
    559             if self.spatial is False:
    560                 raise 'Triangles and vertex_coordinates must be specified'
    561            
     625            #no longer true. sts files have spatial = True but
     626            #if self.spatial is False:
     627            #    raise 'Triangles and vertex_coordinates must be specified'
     628            #
    562629            try:
    563630                self.interpolation_points = interpolation_points = ensure_numeric(interpolation_points)
    … …  
    569636                raise msg
    570637
    571 
    572             # Check that all interpolation points fall within
    573             # mesh boundary as defined by triangles and vertex_coordinates.
    574             from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
    575             from anuga.utilities.polygon import outside_polygon           
    576 
    577             # Create temporary mesh object from mesh info passed
    578             # into this function.
    579             mesh = Mesh(vertex_coordinates, triangles)
    580             mesh_boundary_polygon = mesh.get_boundary_polygon()
    581 
    582            
    583             indices = outside_polygon(interpolation_points,
    584                                       mesh_boundary_polygon)
    585 
    586             # Record result
    587             #self.mesh_boundary_polygon = mesh_boundary_polygon
    588             self.indices_outside_mesh = indices
    589 
    590             # Report
    591             if len(indices) > 0:
    592                 msg = 'Interpolation points in Interpolation function fall '
    593                 msg += 'outside specified mesh. '
    594                 msg += 'Offending points:\n'
    595                 out_interp_pts = []
    596                 for i in indices:
    597                     msg += '%d: %s\n' %(i, interpolation_points[i])
    598                     out_interp_pts.append(ensure_numeric(interpolation_points[i]))
    599 
    600                 if verbose is True:
    601                     import sys
    602                     if sys.platform == 'win32':
    603                         from anuga.utilities.polygon import plot_polygons
    604                         #out_interp_pts = take(interpolation_points,[indices])
    605                         title = 'Interpolation points fall outside specified mesh'
    606                         plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts],
    607                                       ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose)
    608 
    609                 # Joaquim Luis suggested this as an Exception, so
    610                 # that the user can now what the problem is rather than
    611                 # looking for NaN's. However, NANs are handy as they can
    612                 # be ignored leaving good points for continued processing.
    613                 if verbose:
    614                     print msg
    615                 #raise Exception(msg)
     638            if triangles is not None and vertex_coordinates is not None:
     639                # Check that all interpolation points fall within
     640                # mesh boundary as defined by triangles and vertex_coordinates.
     641                from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
     642                from anuga.utilities.polygon import outside_polygon           
     643
     644                # Create temporary mesh object from mesh info passed
     645                # into this function.
     646                mesh = Mesh(vertex_coordinates, triangles)
     647                mesh_boundary_polygon = mesh.get_boundary_polygon()
     648
     649           
     650                indices = outside_polygon(interpolation_points,
     651                                          mesh_boundary_polygon)
     652
     653                # Record result
     654                #self.mesh_boundary_polygon = mesh_boundary_polygon
     655                self.indices_outside_mesh = indices
     656
     657                # Report
     658                if len(indices) > 0:
     659                    msg = 'Interpolation points in Interpolation function fall '
     660                    msg += 'outside specified mesh. '
     661                    msg += 'Offending points:\n'
     662                    out_interp_pts = []
     663                    for i in indices:
     664                        msg += '%d: %s\n' %(i, interpolation_points[i])
     665                        out_interp_pts.append(ensure_numeric(interpolation_points[i]))
     666
     667                    if verbose is True:
     668                        import sys
     669                        if sys.platform == 'win32':
     670                            from anuga.utilities.polygon import plot_polygons
     671                            #out_interp_pts = take(interpolation_points,[indices])
     672                            title = 'Interpolation points fall outside specified mesh'
     673                            plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts],
     674                                          ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose)
     675
     676                    # Joaquim Luis suggested this as an Exception, so
     677                    # that the user can now what the problem is rather than
     678                    # looking for NaN's. However, NANs are handy as they can
     679                    # be ignored leaving good points for continued processing.
     680                    if verbose:
     681                        print msg
     682                    #raise Exception(msg)
     683            elif triangles is None and vertex_coordinates is not None:#jj
     684                #Dealing with sts file
     685                pass
     686            else:
     687                msg = 'Sww file function requires both triangles and vertex_coordinates. sts file file function requires the later.'
     688                raise Exception(msg)
    616689
    617690            # Plot boundary and interpolation points
    … …  
    632705            # Build interpolator
    633706            if verbose:
    634                 msg = 'Building interpolation matrix from source mesh '
    635                 msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
    636                                                        triangles.shape[0])
     707                if triangles is not None and vertex_coordinates is not None:
     708                    msg = 'Building interpolation matrix from source mesh '
     709                    msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
     710                                                           triangles.shape[0])
     711                elif triangles is None and vertex_coordinates is not None:
     712                    msg = 'Building interpolation matrix from source points'
     713               
    637714                print msg
    638715
    … …  
    665742                        print '    quantity %s, size=%d' %(name, len(Q))
    666743                       
    667                     # Interpolate   
    668                     result = interpol.interpolate(Q,
    669                                                   point_coordinates=\
    670                                                   self.interpolation_points,
    671                                                   verbose=False) # Don't clutter
     744                    # Interpolate
     745                    if triangles is not None and vertex_coordinates is not None:   
     746                        result = interpol.interpolate(Q,
     747                                                      point_coordinates=\
     748                                                      self.interpolation_points,
     749                                                      verbose=False) # Don't clutter
     750                    elif triangles is None and vertex_coordinates is not None:
     751                        result=interpol.interpolate_polyline(Q,vertex_coordinates,point_coordinates=self.interpolation_points)
    672752
    673753                    #assert len(result), len(interpolation_points)

  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5358 r5361  
    62896289
    62906290    def test_urs2sts(self):
    6291         """tide = 1
    6292         time_step_count = 3
    6293         time_step = 2
    6294         lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
    6295         depth=20
    6296         ha=2
    6297         ua=5
    6298         va=-10 #-ve added to take into account mux file format where south
    6299                # is positive.
    6300         """
    6301         #tide = 1
    63026291        tide=0
    63036292        time_step_count = 3
    … …  
    64076396        self.delete_mux(files)
    64086397        os.remove(sts_file)
     6398
     6399    def test_file_boundary_sts(self):
     6400        from anuga.shallow_water import Domain
     6401        from anuga.shallow_water import Reflective_boundary
     6402        from anuga.shallow_water import Dirichlet_boundary
     6403        from anuga.shallow_water import File_boundary
     6404        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     6405        from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
     6406        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
     6407        tide=0.
     6408        time_step_count = 5
     6409        time_step = 2
     6410        lat_long_points =bounding_polygon[0:3]
     6411        n=len(lat_long_points)
     6412        first_tstep=ones(n,Int)
     6413        last_tstep=(time_step_count)*ones(n,Int)
     6414        gauge_depth=20*ones(n,Float)
     6415        ha=2*ones((n,time_step_count),Float)
     6416        ua=10*ones((n,time_step_count),Float)
     6417        va=-10*ones((n,time_step_count),Float)
     6418        base_name, files = self.write_mux2(lat_long_points,
     6419                                   time_step_count, time_step,
     6420                                   first_tstep, last_tstep,
     6421                                   depth=gauge_depth,
     6422                                   ha=ha,
     6423                                   ua=ua,
     6424                                   va=va)
     6425
     6426        sts_file=base_name
     6427        urs2sts(base_name,sts_file,mean_stage=tide,verbose=False)
     6428        self.delete_mux(files)
     6429
     6430        #print 'start create mesh from regions'
     6431        for i in range(len(bounding_polygon)):
     6432            zone,bounding_polygon[i][0],bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],bounding_polygon[i][1])
     6433        extent_res=1000000
     6434        meshname = 'urs_test_mesh' + '.tsh'
     6435        interior_regions=None
     6436        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
     6437        create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags,
     6438                         maximum_triangle_area=extent_res,filename=meshname,
     6439                         interior_regions=interior_regions,verbose=False)
     6440       
     6441        domain_fbound = pmesh_to_domain_instance(meshname, Domain)
     6442        domain_fbound.set_quantity('stage', tide)
     6443        Bf = File_boundary(sts_file+'.sts', domain_fbound)
     6444        Br = Reflective_boundary(domain_fbound)
     6445        Bd=Dirichlet_boundary([2.0,220,-220])
     6446        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
     6447        finaltime=time_step*(time_step_count-1)
     6448        yieldstep=time_step
     6449        temp_fbound=zeros(int(finaltime/yieldstep)+1,Float)
     6450        i=0
     6451        for t in domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,
     6452                                      skip_initial_step = False):
     6453            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
     6454            i+=1
     6455       
     6456        domain_drchlt = pmesh_to_domain_instance(meshname, Domain)
     6457        domain_drchlt.set_quantity('stage', tide)
     6458        Br = Reflective_boundary(domain_drchlt)
     6459        Bd=Dirichlet_boundary([2.0,220,-220])
     6460        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
     6461        temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float)
     6462        i=0
     6463        for t in domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,
     6464                                      skip_initial_step = False):
     6465            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
     6466            i+=1
     6467
     6468        assert temp_fbound-temp_drchlt<epsilon
     6469        os.remove(sts_file+'.sts')
     6470        os.remove(meshname)
    64096471
    64106472    def test_lon_lat2grid(self):
Note: See TracChangeset for help on using the changeset viewer.