Changeset 1884 for inundation/pyvolution


Ignore:
Timestamp:
Oct 10, 2005, 11:41:26 AM (18 years ago)
Author:
ole
Message:

Made interpolation_points in file_function absolute UTM coordinates and wrote test

Location:
inundation/pyvolution
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/least_squares.py

    r1882 r1884  
    939939        from Numeric import array, zeros, Float, alltrue, concatenate,\
    940940             reshape, ArrayType
     941
    941942
    942943        from util import mean, ensure_numeric
  • inundation/pyvolution/test_util.py

    r1854 r1884  
    472472        os.remove(filename)
    473473
    474        
     474
     475
     476    def test_spatio_temporal_file_function_different_origin(self):
     477        """Test that spatio temporal file function performs the correct
     478        interpolations in both time and space where space is offset by
     479        xllcorner and yllcorner
     480        NetCDF version (x,y,t dependency)       
     481        """
     482        import time
     483
     484        #Create sww file of simple propagation from left to right
     485        #through rectangular domain
     486        from shallow_water import Domain, Dirichlet_boundary
     487        from mesh_factory import rectangular
     488        from Numeric import take, concatenate, reshape
     489
     490
     491        from coordinate_transforms.geo_reference import Geo_reference
     492        xllcorner = 2048
     493        yllcorner = 11000
     494        zone = 2
     495
     496        #Create basic mesh and shallow water domain
     497        points, vertices, boundary = rectangular(3, 3)
     498        domain1 = Domain(points, vertices, boundary,
     499                         geo_reference = Geo_reference(xllcorner = xllcorner,
     500                                                       yllcorner = yllcorner))
     501       
     502
     503        from util import mean
     504        domain1.reduction = mean
     505        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
     506                              # only one value.
     507
     508        domain1.default_order = 2
     509        domain1.store = True
     510        domain1.set_datadir('.')
     511        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
     512        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
     513
     514        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
     515        domain1.set_quantity('elevation', 0)
     516        domain1.set_quantity('friction', 0)
     517        domain1.set_quantity('stage', 0)
     518
     519        # Boundary conditions
     520        B0 = Dirichlet_boundary([0,0,0])
     521        B6 = Dirichlet_boundary([0.6,0,0])
     522        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
     523        domain1.check_integrity()
     524
     525        finaltime = 8
     526        #Evolution
     527        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
     528            pass
     529            #domain1.write_time()
     530
     531
     532        #Now read data from sww and check
     533        from Scientific.IO.NetCDF import NetCDFFile
     534        filename = domain1.get_name() + '.' + domain1.format
     535        fid = NetCDFFile(filename)
     536
     537        x = fid.variables['x'][:]
     538        y = fid.variables['y'][:]
     539        stage = fid.variables['stage'][:]
     540        xmomentum = fid.variables['xmomentum'][:]
     541        ymomentum = fid.variables['ymomentum'][:]
     542        time = fid.variables['time'][:]
     543
     544        #Take stage vertex values at last timestep on diagonal
     545        #Diagonal is identified by vertices: 0, 5, 10, 15
     546
     547        timestep = len(time)-1 #Last timestep
     548        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     549        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     550        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     551        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     552
     553        #Reference interpolated values at midpoints on diagonal at
     554        #this timestep are
     555        r0 = (D[0] + D[1])/2
     556        r1 = (D[1] + D[2])/2
     557        r2 = (D[2] + D[3])/2
     558
     559        #And the midpoints are found now
     560        Dx = take(reshape(x, (16,1)), [0,5,10,15])
     561        Dy = take(reshape(y, (16,1)), [0,5,10,15])
     562
     563        diag = concatenate( (Dx, Dy), axis=1)
     564        d_midpoints = (diag[1:] + diag[:-1])/2
     565
     566
     567        #Adjust for georef - make interpolation points absolute
     568        d_midpoints[:,0] += xllcorner
     569        d_midpoints[:,1] += yllcorner               
     570
     571        #Let us see if the file function can find the correct
     572        #values at the midpoints at the last timestep:
     573        f = file_function(filename, domain1,
     574                          interpolation_points = d_midpoints)
     575
     576        q = f(timestep/10., point_id=0); assert allclose(r0, q)
     577        q = f(timestep/10., point_id=1); assert allclose(r1, q)
     578        q = f(timestep/10., point_id=2); assert allclose(r2, q)
     579
     580
     581        ##################
     582        #Now do the same for the first timestep
     583
     584        timestep = 0 #First timestep
     585        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     586        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     587        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     588        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     589
     590        #Reference interpolated values at midpoints on diagonal at
     591        #this timestep are
     592        r0 = (D[0] + D[1])/2
     593        r1 = (D[1] + D[2])/2
     594        r2 = (D[2] + D[3])/2
     595
     596        #Let us see if the file function can find the correct
     597        #values
     598        q = f(0, point_id=0); assert allclose(r0, q)
     599        q = f(0, point_id=1); assert allclose(r1, q)
     600        q = f(0, point_id=2); assert allclose(r2, q)
     601
     602
     603        ##################
     604        #Now do it again for a timestep in the middle
     605
     606        timestep = 33
     607        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     608        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     609        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     610        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     611
     612        #Reference interpolated values at midpoints on diagonal at
     613        #this timestep are
     614        r0 = (D[0] + D[1])/2
     615        r1 = (D[1] + D[2])/2
     616        r2 = (D[2] + D[3])/2
     617
     618        q = f(timestep/10., point_id=0); assert allclose(r0, q)
     619        q = f(timestep/10., point_id=1); assert allclose(r1, q)
     620        q = f(timestep/10., point_id=2); assert allclose(r2, q)
     621
     622
     623        ##################
     624        #Now check temporal interpolation
     625        #Halfway between timestep 15 and 16
     626
     627        timestep = 15
     628        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     629        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     630        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     631        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     632
     633        #Reference interpolated values at midpoints on diagonal at
     634        #this timestep are
     635        r0_0 = (D[0] + D[1])/2
     636        r1_0 = (D[1] + D[2])/2
     637        r2_0 = (D[2] + D[3])/2
     638
     639        #
     640        timestep = 16
     641        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
     642        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
     643        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     644        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
     645
     646        #Reference interpolated values at midpoints on diagonal at
     647        #this timestep are
     648        r0_1 = (D[0] + D[1])/2
     649        r1_1 = (D[1] + D[2])/2
     650        r2_1 = (D[2] + D[3])/2
     651
     652        # The reference values are
     653        r0 = (r0_0 + r0_1)/2
     654        r1 = (r1_0 + r1_1)/2
     655        r2 = (r2_0 + r2_1)/2
     656
     657        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
     658        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
     659        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
     660
     661        ##################
     662        #Finally check interpolation 2 thirds of the way
     663        #between timestep 15 and 16
     664
     665        # The reference values are
     666        r0 = (r0_0 + 2*r0_1)/3
     667        r1 = (r1_0 + 2*r1_1)/3
     668        r2 = (r2_0 + 2*r2_1)/3
     669
     670        #And the file function gives
     671        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
     672        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
     673        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
     674
     675        fid.close()
     676        import os
     677        os.remove(filename)
     678
     679       
     680
    475681
    476682    def test_spatio_temporal_file_function_time(self):
  • inundation/pyvolution/util.py

    r1859 r1884  
    155155    All times are assumed to be in UTC
    156156
    157     All spatial information is assumed to be in UTM coordinates.
     157    All spatial information is assumed to be in absolute UTM coordinates.
    158158
    159159    See Interpolation function for further documentation
     
    170170    #- domain's georef
    171171    #- sww file's georef
    172     #- interpolation points georef
     172    #- interpolation points as absolute UTM coordinates
    173173
    174174
     
    228228    from config import time_format
    229229    from Scientific.IO.NetCDF import NetCDFFile
    230     from Numeric import array, zeros, Float, alltrue, concatenate, reshape       
     230    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
     231    from util import ensure_numeric   
    231232
    232233    #Open NetCDF file
     
    251252
    252253
     254    if interpolation_points is not None:
     255        interpolation_points = ensure_numeric(interpolation_points)
     256
     257
     258
    253259    #Now assert that requested quantitites (and the independent ones)
    254260    #are present in file
     
    302308        y = reshape(y, (len(y),1))
    303309        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
     310
     311        if interpolation_points is not None:
     312            #Adjust for georef
     313            interpolation_points[:,0] -= xllcorner
     314            interpolation_points[:,1] -= yllcorner       
     315       
     316
    304317
    305318
Note: See TracChangeset for help on using the changeset viewer.