Python Domain.geo_reference方法代码示例

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import geo_reference [as 别名]
    def test_sww2domain1(self):
        ################################################
        #Create a test domain, and evolve and save it.
        ################################################
        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        #Create basic mesh

        yiel=0.01
        points, vertices, boundary = rectangular(10,10)

        #print "=============== boundary rect ======================="
        #print boundary

        #Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.geo_reference = Geo_reference(56,11,11)
        domain.smooth = False
        domain.store = True
        domain.set_name('bedslope')
        domain.default_order=2
        #Bed-slope and friction
        domain.set_quantity('elevation', lambda x,y: -x/3)
        domain.set_quantity('friction', 0.1)
        # Boundary conditions
        from math import sin, pi
        Br = Reflective_boundary(domain)
        Bt = Transmissive_boundary(domain)
        Bd = Dirichlet_boundary([0.2,0.,0.])
        Bw = Time_boundary(domain=domain,function=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])

        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})

        domain.quantities_to_be_stored['xmomentum'] = 2
        domain.quantities_to_be_stored['ymomentum'] = 2
        #Initial condition
        h = 0.05
        elevation = domain.quantities['elevation'].vertex_values
        domain.set_quantity('stage', elevation + h)

        domain.check_integrity()
        #Evolution
        #domain.tight_slope_limiters = 1
        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
            #domain.write_time()
            pass

        #print boundary


        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
        domain2 = load_sww_as_domain(filename, None, fail_if_NaN=False,
                                        verbose=self.verbose)

        # Unfortunately we loss the boundaries top, bottom, left and right,
        # they are now all lumped into "exterior"

        #print "=============== boundary domain2 ======================="
        #print domain2.boundary
        

        #print domain2.get_boundary_tags()
        
        #points, vertices, boundary = rectangular(15,15)
        #domain2.boundary = boundary
        ###################
        ##NOW TEST IT!!!
        ###################

        os.remove(filename)

        bits = ['vertex_coordinates']
        for quantity in ['stage']:
            bits.append('get_quantity("%s").get_integral()' % quantity)
            bits.append('get_quantity("%s").get_values()' % quantity)

        for bit in bits:
            #print 'testing that domain.'+bit+' has been restored'
            #print bit
            #print 'done'
            #print eval('domain.'+bit)
            #print eval('domain2.'+bit)
            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))

        ######################################
        #Now evolve them both, just to be sure
        ######################################x
        from time import sleep

        final = .1
        domain.set_quantity('friction', 0.1)
        domain.store = False
        domain.set_boundary({'exterior': Bd, 'left' : Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})


        for t in domain.evolve(yieldstep = yiel, finaltime = final):
            #domain.write_time()
            pass

#.........这里部分代码省略.........

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import geo_reference [as 别名]
    def test_sww2pts_centroids_de0(self):
        """Test that sww information can be converted correctly to pts data at specified coordinates
        - in this case, the centroids.
        """

        import time, os
        from anuga.file.netcdf import NetCDFFile
        # Used for points that lie outside mesh
        NODATA_value = 1758323

        # Setup
        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create shallow water domain
        domain = Domain(*rectangular(2, 2))

        B = Transmissive_boundary(domain)
        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})

        domain.set_name('datatest_de0')

        ptsfile = domain.get_name() + '_elevation.pts'
        swwfile = domain.get_name() + '.sww'

        domain.set_datadir('.')
        domain.format = 'sww'
        domain.set_quantity('elevation', lambda x,y: -x-y)

        domain.geo_reference = Geo_reference(56,308500,6189000)

        sww = SWW_file(domain)
        sww.store_connectivity()
        sww.store_timestep()

        #self.domain.tight_slope_limiters = 1
        domain.evolve_to_end(finaltime = 0.01)
        sww.store_timestep()

        # Check contents in NetCDF
        fid = NetCDFFile(sww.filename, netcdf_mode_r)

        # Get the variables
        x = fid.variables['x'][:]
        y = fid.variables['y'][:]
        elevation = fid.variables['elevation'][:]
        time = fid.variables['time'][:]
        stage = fid.variables['stage'][:]

        volumes = fid.variables['volumes'][:]


        # Invoke interpolation for vertex points       
        points = num.concatenate( (x[:,num.newaxis],y[:,num.newaxis]), axis=1 )
        points = num.ascontiguousarray(points)
        sww2pts(domain.get_name() + '.sww',
                quantity = 'elevation',
                data_points = points,
                NODATA_value = NODATA_value)
        ref_point_values = elevation
        point_values = Geospatial_data(ptsfile).get_attributes()
        #print 'P', point_values
        #print 'Ref', ref_point_values        
        assert num.allclose(point_values, ref_point_values)        



        # Invoke interpolation for centroids
        points = domain.get_centroid_coordinates()
        #print points
        sww2pts(domain.get_name() + '.sww',
                quantity = 'elevation',
                data_points = points,
                NODATA_value = NODATA_value)
        #ref_point_values = [-0.5, -0.5, -1, -1, -1, -1, -1.5, -1.5]   #At centroids

        ref_point_values = [-0.77777777, -0.77777777, -0.99999998, -0.99999998, 
                             -0.99999998, -0.99999998, -1.22222221, -1.22222221]
        point_values = Geospatial_data(ptsfile).get_attributes()
        #print 'P', point_values
        #print 'Ref', ref_point_values        
        assert num.allclose(point_values, ref_point_values)        

        fid.close()

        #Cleanup
        os.remove(sww.filename)
        os.remove(ptsfile)

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import geo_reference [as 别名]

#.........这里部分代码省略.........
        location = get_maximum_inundation_location(swwfile)
        #print 'Runup, location', runup, location
        assert num.allclose(runup, 4.66666666667)
        assert num.allclose(location[0], 46.666668) 
               
        # Check final runup
        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
        #print 'Runup, location:',runup, location

        assert num.allclose(runup, 3.81481488546)
        assert num.allclose(location[0], 51.666668)

        # Check runup restricted to a polygon
        p = [[50,1], [99,1], [99,49], [50,49]]
        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
        location = get_maximum_inundation_location(swwfile, polygon=p)
        #print runup, location

        assert num.allclose(runup, 3.81481488546) 
        assert num.allclose(location[0], 51.6666666)                

        # Check that mimimum_storable_height works
        fid = NetCDFFile(swwfile, netcdf_mode_r) # Open existing file
        
        stage = fid.variables['stage_c'][:]
        z = fid.variables['elevation_c'][:]
        xmomentum = fid.variables['xmomentum_c'][:]
        ymomentum = fid.variables['ymomentum_c'][:]
        
        for i in range(stage.shape[0]):
            h = stage[i]-z # depth vector at time step i
            
            # Check every node location
            for j in range(stage.shape[1]):
                # Depth being either exactly zero implies
                # momentum being zero.
                # Or else depth must be greater than or equal to
                # the minimal storable height
                if h[j] == 0.0:
                    assert xmomentum[i,j] == 0.0
                    assert ymomentum[i,j] == 0.0                
                else:
                    assert h[j] >= 0.0
        
        fid.close()

        # Cleanup
        os.remove(swwfile)
        


        #------------- Now the same with georeferencing

        domain.time=0.0
        E = 308500
        N = 6189000
        #E = N = 0
        domain.geo_reference = Geo_reference(56, E, N)

        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
        domain.set_quantity('stage', -6)
        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})

        for t in domain.evolve(yieldstep=1, finaltime = 50):
            pass

        # Check maximal runup
        runup = get_maximum_inundation_elevation(swwfile)
        location = get_maximum_inundation_location(swwfile)

        #print runup, location

        assert num.allclose(runup,4.66666666667)
        assert num.allclose(location[0], 308546.66) 

        # Check final runup
        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
        
        #print runup, location
        #1.66666666667 [308561.66, 6189006.5]

        assert num.allclose(runup, 3.81481488546)
        assert num.allclose(location[0], 308551.66)

        # Check runup restricted to a polygon
        p = num.array([[50,1], [99,1], [99,49], [50,49]], num.int) + num.array([E, N], num.int)      #array default#

        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
        location = get_maximum_inundation_location(swwfile, polygon=p)

        #print runup, location

        assert num.allclose(runup, 3.81481488546)
        assert num.allclose(location[0], 308551.66)                


        # Cleanup
        os.remove(swwfile)