#!/usr/bin/env python import unittest from Numeric import zeros, array, allclose from math import sqrt, pi from util import * from config import epsilon from data_manager import timefile2netcdf def test_function(x, y): return x+y class Test_Util(unittest.TestCase): def setUp(self): pass def tearDown(self): pass def test_gradient(self): x0 = 0.0; y0 = 0.0; z0 = 0.0 x1 = 1.0; y1 = 0.0; z1 = -1.0 x2 = 0.0; y2 = 1.0; z2 = 0.0 zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) assert zx == -1.0 assert zy == 0.0 def test_gradient_more(self): x0 = 2.0/3; y0 = 2.0/3 x1= 8.0/3; y1 = 2.0/3 x2 = 2.0/3; y2 = 8.0/3 q0 = 2.0+2.0/3 q1 = 8.0+2.0/3 q2 = 2.0+8.0/3 #Gradient of fitted pwl surface a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) assert abs(a - 3.0) < epsilon assert abs(b - 1.0) < epsilon def test_gradient2(self): """Test two-point gradient """ x0 = 5.0; y0 = 5.0; z0 = 10.0 x1 = 8.0; y1 = 2.0; z1 = 1.0 x2 = 8.0; y2 = 8.0; z2 = 10.0 #Reference zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) a, b = gradient2(x0, y0, x1, y1, z0, z1) assert zx == a assert zy == b z2_computed = z0 + a*(x2-x0) + b*(y2-y0) assert z2_computed == z2 def test_gradient2_more(self): """Test two-point gradient more """ x0 = 2.0; y0 = 2.0 x1 = 8.0; y1 = 3.0 x2 = 1.0; y2 = 8.0 q0 = 2.0 q1 = 8.0 q2 = q0 #Gradient of fitted pwl surface a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) a, b = gradient2(x0, y0, x1, y1, q0, q1) assert a == a_ref assert b == b_ref def test_that_C_extension_compiles(self): FN = 'util_ext.c' try: import util_ext except: from compile import compile try: compile(FN) except: raise 'Could not compile %s' %FN else: import util_ext def test_gradient_C_extension(self): from util_ext import gradient as gradient_c x0 = 2.0/3; y0 = 2.0/3 x1= 8.0/3; y1 = 2.0/3 x2 = 2.0/3; y2 = 8.0/3 q0 = 2.0+2.0/3 q1 = 8.0+2.0/3 q2 = 2.0+8.0/3 #Gradient of fitted pwl surface a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2) assert abs(a - 3.0) < epsilon assert abs(b - 1.0) < epsilon def test_gradient_C_extension3(self): from util_ext import gradient as gradient_c from RandomArray import uniform, seed seed(17, 53) x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6) q0 = uniform(0.0, 10.0, 4) q1 = uniform(1.0, 3.0, 4) q2 = uniform(7.0, 20.0, 4) for i in range(4): #Gradient of fitted pwl surface from util import gradient_python a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0[i], q1[i], q2[i]) #print a_ref, b_ref a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0[i], q1[i], q2[i]) #print a, a_ref, b, b_ref assert abs(a - a_ref) < epsilon assert abs(b - b_ref) < epsilon #Geometric #def test_distance(self): # from util import distance# # # self.failUnless( distance([4,2],[7,6]) == 5.0, # 'Distance is wrong!') # self.failUnless( allclose(distance([7,6],[9,8]), 2.82842712475), # 'distance is wrong!') # self.failUnless( allclose(distance([9,8],[4,2]), 7.81024967591), # 'distance is wrong!') # # self.failUnless( distance([9,8],[4,2]) == distance([4,2],[9,8]), # 'distance is wrong!') def test_angle(self): from util import angle assert allclose(angle([1.0, 1.0])/pi*180, 45.0) def test_anglediff(self): from util import anglediff assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0) def test_ensure_numeric(self): from util import ensure_numeric from Numeric import ArrayType, Float, array A = [1,2,3,4] B = ensure_numeric(A) assert type(B) == ArrayType assert B.typecode() == 'l' assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4 A = [1,2,3.14,4] B = ensure_numeric(A) assert type(B) == ArrayType assert B.typecode() == 'd' assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4 A = [1,2,3,4] B = ensure_numeric(A, Float) assert type(B) == ArrayType assert B.typecode() == 'd' assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0 A = [1,2,3,4] B = ensure_numeric(A, Float) assert type(B) == ArrayType assert B.typecode() == 'd' assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0 A = array([1,2,3,4]) B = ensure_numeric(A) assert type(B) == ArrayType assert B.typecode() == 'l' assert A == B assert A is B #Same object A = array([1,2,3,4]) B = ensure_numeric(A, Float) assert type(B) == ArrayType assert B.typecode() == 'd' assert A == B assert A is not B #Not the same object def test_file_function_time1(self): """Test that File function interpolates correctly between given times. No x,y dependency here. """ #Write file import os, time from config import time_format from math import sin, pi #Typical ASCII file finaltime = 1200 filename = 'test_file_function' fid = open(filename + '.txt', 'w') start = time.mktime(time.strptime('2000', '%Y')) dt = 60 #One minute intervals t = 0.0 while t <= finaltime: t_string = time.strftime(time_format, time.gmtime(t+start)) fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600))) t += dt fid.close() #Convert ASCII file to NetCDF (Which is what we really like!) timefile2netcdf(filename) #Create file function from time series F = file_function(filename + '.tms', quantities = ['Attribute0', 'Attribute1', 'Attribute2']) #Now try interpolation for i in range(20): t = i*10 q = F(t) #Exact linear intpolation assert allclose(q[0], 2*t) if i%6 == 0: assert allclose(q[1], t**2) assert allclose(q[2], sin(t*pi/600)) #Check non-exact t = 90 #Halfway between 60 and 120 q = F(t) assert allclose( (120**2 + 60**2)/2, q[1] ) assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] ) t = 100 #Two thirds of the way between between 60 and 120 q = F(t) assert allclose( 2*120**2/3 + 60**2/3, q[1] ) assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] ) os.remove(filename + '.txt') os.remove(filename + '.tms') def test_spatio_temporal_file_function(self): """Test that spatio temporal file function performs the correct interpolations in both time and space NetCDF version (x,y,t dependency) """ import time #Create sww file of simple propagation from left to right #through rectangular domain from shallow_water import Domain, Dirichlet_boundary from mesh_factory import rectangular from Numeric import take, concatenate, reshape #Create basic mesh and shallow water domain points, vertices, boundary = rectangular(3, 3) domain1 = Domain(points, vertices, boundary) from util import mean domain1.reduction = mean domain1.smooth = True #NOTE: Mimic sww output where each vertex has # only one value. domain1.default_order = 2 domain1.store = True domain1.set_datadir('.') domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self))) domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] #Bed-slope, friction and IC at vertices (and interpolated elsewhere) domain1.set_quantity('elevation', 0) domain1.set_quantity('friction', 0) domain1.set_quantity('stage', 0) # Boundary conditions B0 = Dirichlet_boundary([0,0,0]) B6 = Dirichlet_boundary([0.6,0,0]) domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0}) domain1.check_integrity() finaltime = 8 #Evolution for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime): pass #domain1.write_time() #Now read data from sww and check from Scientific.IO.NetCDF import NetCDFFile filename = domain1.get_name() + '.' + domain1.format fid = NetCDFFile(filename) x = fid.variables['x'][:] y = fid.variables['y'][:] stage = fid.variables['stage'][:] xmomentum = fid.variables['xmomentum'][:] ymomentum = fid.variables['ymomentum'][:] time = fid.variables['time'][:] #Take stage vertex values at last timestep on diagonal #Diagonal is identified by vertices: 0, 5, 10, 15 timestep = len(time)-1 #Last timestep d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0 = (D[0] + D[1])/2 r1 = (D[1] + D[2])/2 r2 = (D[2] + D[3])/2 #And the midpoints are found now Dx = take(reshape(x, (16,1)), [0,5,10,15]) Dy = take(reshape(y, (16,1)), [0,5,10,15]) diag = concatenate( (Dx, Dy), axis=1) d_midpoints = (diag[1:] + diag[:-1])/2 #Let us see if the file function can find the correct #values at the midpoints at the last timestep: f = file_function(filename, domain1, interpolation_points = d_midpoints) q = f(timestep/10., point_id=0); assert allclose(r0, q) q = f(timestep/10., point_id=1); assert allclose(r1, q) q = f(timestep/10., point_id=2); assert allclose(r2, q) ################## #Now do the same for the first timestep timestep = 0 #First timestep d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0 = (D[0] + D[1])/2 r1 = (D[1] + D[2])/2 r2 = (D[2] + D[3])/2 #Let us see if the file function can find the correct #values q = f(0, point_id=0); assert allclose(r0, q) q = f(0, point_id=1); assert allclose(r1, q) q = f(0, point_id=2); assert allclose(r2, q) ################## #Now do it again for a timestep in the middle timestep = 33 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0 = (D[0] + D[1])/2 r1 = (D[1] + D[2])/2 r2 = (D[2] + D[3])/2 q = f(timestep/10., point_id=0); assert allclose(r0, q) q = f(timestep/10., point_id=1); assert allclose(r1, q) q = f(timestep/10., point_id=2); assert allclose(r2, q) ################## #Now check temporal interpolation #Halfway between timestep 15 and 16 timestep = 15 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0_0 = (D[0] + D[1])/2 r1_0 = (D[1] + D[2])/2 r2_0 = (D[2] + D[3])/2 # timestep = 16 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0_1 = (D[0] + D[1])/2 r1_1 = (D[1] + D[2])/2 r2_1 = (D[2] + D[3])/2 # The reference values are r0 = (r0_0 + r0_1)/2 r1 = (r1_0 + r1_1)/2 r2 = (r2_0 + r2_1)/2 q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q) q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q) q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q) ################## #Finally check interpolation 2 thirds of the way #between timestep 15 and 16 # The reference values are r0 = (r0_0 + 2*r0_1)/3 r1 = (r1_0 + 2*r1_1)/3 r2 = (r2_0 + 2*r2_1)/3 #And the file function gives q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q) q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q) q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q) fid.close() import os os.remove(filename) def test_spatio_temporal_file_function_different_origin(self): """Test that spatio temporal file function performs the correct interpolations in both time and space where space is offset by xllcorner and yllcorner NetCDF version (x,y,t dependency) """ import time #Create sww file of simple propagation from left to right #through rectangular domain from shallow_water import Domain, Dirichlet_boundary from mesh_factory import rectangular from Numeric import take, concatenate, reshape from coordinate_transforms.geo_reference import Geo_reference xllcorner = 2048 yllcorner = 11000 zone = 2 #Create basic mesh and shallow water domain points, vertices, boundary = rectangular(3, 3) domain1 = Domain(points, vertices, boundary, geo_reference = Geo_reference(xllcorner = xllcorner, yllcorner = yllcorner)) from util import mean domain1.reduction = mean domain1.smooth = True #NOTE: Mimic sww output where each vertex has # only one value. domain1.default_order = 2 domain1.store = True domain1.set_datadir('.') domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self))) domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] #Bed-slope, friction and IC at vertices (and interpolated elsewhere) domain1.set_quantity('elevation', 0) domain1.set_quantity('friction', 0) domain1.set_quantity('stage', 0) # Boundary conditions B0 = Dirichlet_boundary([0,0,0]) B6 = Dirichlet_boundary([0.6,0,0]) domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0}) domain1.check_integrity() finaltime = 8 #Evolution for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime): pass #domain1.write_time() #Now read data from sww and check from Scientific.IO.NetCDF import NetCDFFile filename = domain1.get_name() + '.' + domain1.format fid = NetCDFFile(filename) x = fid.variables['x'][:] y = fid.variables['y'][:] stage = fid.variables['stage'][:] xmomentum = fid.variables['xmomentum'][:] ymomentum = fid.variables['ymomentum'][:] time = fid.variables['time'][:] #Take stage vertex values at last timestep on diagonal #Diagonal is identified by vertices: 0, 5, 10, 15 timestep = len(time)-1 #Last timestep d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0 = (D[0] + D[1])/2 r1 = (D[1] + D[2])/2 r2 = (D[2] + D[3])/2 #And the midpoints are found now Dx = take(reshape(x, (16,1)), [0,5,10,15]) Dy = take(reshape(y, (16,1)), [0,5,10,15]) diag = concatenate( (Dx, Dy), axis=1) d_midpoints = (diag[1:] + diag[:-1])/2 #Adjust for georef - make interpolation points absolute d_midpoints[:,0] += xllcorner d_midpoints[:,1] += yllcorner #Let us see if the file function can find the correct #values at the midpoints at the last timestep: f = file_function(filename, domain1, interpolation_points = d_midpoints) q = f(timestep/10., point_id=0); assert allclose(r0, q) q = f(timestep/10., point_id=1); assert allclose(r1, q) q = f(timestep/10., point_id=2); assert allclose(r2, q) ################## #Now do the same for the first timestep timestep = 0 #First timestep d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0 = (D[0] + D[1])/2 r1 = (D[1] + D[2])/2 r2 = (D[2] + D[3])/2 #Let us see if the file function can find the correct #values q = f(0, point_id=0); assert allclose(r0, q) q = f(0, point_id=1); assert allclose(r1, q) q = f(0, point_id=2); assert allclose(r2, q) ################## #Now do it again for a timestep in the middle timestep = 33 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0 = (D[0] + D[1])/2 r1 = (D[1] + D[2])/2 r2 = (D[2] + D[3])/2 q = f(timestep/10., point_id=0); assert allclose(r0, q) q = f(timestep/10., point_id=1); assert allclose(r1, q) q = f(timestep/10., point_id=2); assert allclose(r2, q) ################## #Now check temporal interpolation #Halfway between timestep 15 and 16 timestep = 15 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0_0 = (D[0] + D[1])/2 r1_0 = (D[1] + D[2])/2 r2_0 = (D[2] + D[3])/2 # timestep = 16 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) D = concatenate( (d_stage, d_uh, d_vh), axis=1) #Reference interpolated values at midpoints on diagonal at #this timestep are r0_1 = (D[0] + D[1])/2 r1_1 = (D[1] + D[2])/2 r2_1 = (D[2] + D[3])/2 # The reference values are r0 = (r0_0 + r0_1)/2 r1 = (r1_0 + r1_1)/2 r2 = (r2_0 + r2_1)/2 q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q) q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q) q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q) ################## #Finally check interpolation 2 thirds of the way #between timestep 15 and 16 # The reference values are r0 = (r0_0 + 2*r0_1)/3 r1 = (r1_0 + 2*r1_1)/3 r2 = (r2_0 + 2*r2_1)/3 #And the file function gives q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q) q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q) q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q) fid.close() import os os.remove(filename) def test_spatio_temporal_file_function_time(self): """Test that File function interpolates correctly between given times. NetCDF version (x,y,t dependency) """ #Create NetCDF (sww) file to be read # x: 0, 5, 10, 15 # y: -20, -10, 0, 10 # t: 0, 60, 120, ...., 1200 # # test quantities (arbitrary but non-trivial expressions): # # stage = 3*x - y**2 + 2*t # xmomentum = exp( -((x-7)**2 + (y+5)**2)/20 ) * t**2 # ymomentum = x**2 + y**2 * sin(t*pi/600) #NOTE: Nice test that may render some of the others redundant. import os, time from config import time_format from Numeric import sin, pi, exp from mesh_factory import rectangular from shallow_water import Domain import data_manager finaltime = 1200 filename = 'test_file_function' #Create a domain to hold test grid #(0:15, -20:10) points, vertices, boundary =\ rectangular(4, 4, 15, 30, origin = (0, -20)) #print 'Number of elements', len(vertices) domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order = 2 domain.set_datadir('.') domain.set_name(filename) domain.store = True domain.format = 'sww' #Native netcdf visualisation format #print points start = time.mktime(time.strptime('2000', '%Y')) domain.starttime = start #Store structure domain.initialise_storage() #Compute artificial time steps and store dt = 60 #One minute intervals t = 0.0 while t <= finaltime: #Compute quantities f1 = lambda x,y: 3*x - y**2 + 2*t + 4 domain.set_quantity('stage', f1) f2 = lambda x,y: x+y+t**2 domain.set_quantity('xmomentum', f2) f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600) domain.set_quantity('ymomentum', f3) #Store and advance time domain.time = t domain.store_timestep(domain.conserved_quantities) t += dt interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14], [10,-12.5]] #Deliberately set domain.starttime to too early domain.starttime = start - 1 #Create file function F = file_function(filename + '.sww', domain, quantities = domain.conserved_quantities, interpolation_points = interpolation_points) #Check that FF updates fixes domain starttime assert allclose(domain.starttime, start) #Check that domain.starttime isn't updated if later domain.starttime = start + 1 F = file_function(filename + '.sww', domain, quantities = domain.conserved_quantities, interpolation_points = interpolation_points) assert allclose(domain.starttime, start+1) domain.starttime = start #Check linear interpolation in time F = file_function(filename + '.sww', domain, quantities = domain.conserved_quantities, interpolation_points = interpolation_points) for id in range(len(interpolation_points)): x = interpolation_points[id][0] y = interpolation_points[id][1] for i in range(20): t = i*10 k = i%6 if k == 0: q0 = F(t, point_id=id) q1 = F(t+60, point_id=id) q = F(t, point_id=id) #print i, k, t, q #print ' ', q0 #print ' ', q1 #print 's', (k*q1 + (6-k)*q0)/6 #print assert allclose(q, (k*q1 + (6-k)*q0)/6) #Another check of linear interpolation in time for id in range(len(interpolation_points)): q60 = F(60, point_id=id) q120 = F(120, point_id=id) t = 90 #Halfway between 60 and 120 q = F(t, point_id=id) assert allclose( (q120+q60)/2, q ) t = 100 #Two thirds of the way between between 60 and 120 q = F(t, point_id=id) assert allclose(q60/3 + 2*q120/3, q) #Check that domain.starttime isn't updated if later than file starttime but earlier #than file end time delta = 23 domain.starttime = start + delta F = file_function(filename + '.sww', domain, quantities = domain.conserved_quantities, interpolation_points = interpolation_points) assert allclose(domain.starttime, start+delta) #Now try interpolation with delta offset for id in range(len(interpolation_points)): x = interpolation_points[id][0] y = interpolation_points[id][1] for i in range(20): t = i*10 k = i%6 if k == 0: q0 = F(t-delta, point_id=id) q1 = F(t+60-delta, point_id=id) q = F(t-delta, point_id=id) assert allclose(q, (k*q1 + (6-k)*q0)/6) os.remove(filename + '.sww') def test_file_function_time_with_domain(self): """Test that File function interpolates correctly between given times. No x,y dependency here. Use domain with starttime """ #Write file import os, time, calendar from config import time_format from math import sin, pi from domain import Domain finaltime = 1200 filename = 'test_file_function' fid = open(filename + '.txt', 'w') start = time.mktime(time.strptime('2000', '%Y')) dt = 60 #One minute intervals t = 0.0 while t <= finaltime: t_string = time.strftime(time_format, time.gmtime(t+start)) fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600))) t += dt fid.close() #Convert ASCII file to NetCDF (Which is what we really like!) timefile2netcdf(filename) a = [0.0, 0.0] b = [4.0, 0.0] c = [0.0, 3.0] points = [a, b, c] vertices = [[0,1,2]] domain = Domain(points, vertices) #Check that domain.starttime is updated if non-existing F = file_function(filename + '.tms', domain) assert allclose(domain.starttime, start) #Check that domain.starttime is updated if too early domain.starttime = start - 1 F = file_function(filename + '.tms', domain) assert allclose(domain.starttime, start) #Check that domain.starttime isn't updated if later domain.starttime = start + 1 F = file_function(filename + '.tms', domain) assert allclose(domain.starttime, start+1) domain.starttime = start F = file_function(filename + '.tms', domain, quantities = ['Attribute0', 'Attribute1', 'Attribute2']) #print F.T #print F.precomputed_values #print 'F(60)', F(60) #Now try interpolation for i in range(20): t = i*10 q = F(t) #Exact linear intpolation assert allclose(q[0], 2*t) if i%6 == 0: assert allclose(q[1], t**2) assert allclose(q[2], sin(t*pi/600)) #Check non-exact t = 90 #Halfway between 60 and 120 q = F(t) assert allclose( (120**2 + 60**2)/2, q[1] ) assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] ) t = 100 #Two thirds of the way between between 60 and 120 q = F(t) assert allclose( 2*120**2/3 + 60**2/3, q[1] ) assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] ) os.remove(filename + '.tms') os.remove(filename + '.txt') def test_file_function_time_with_domain_different_start(self): """Test that File function interpolates correctly between given times. No x,y dependency here. Use domain with a starttime later than that of file ASCII version """ #Write file import os, time, calendar from config import time_format from math import sin, pi from domain import Domain finaltime = 1200 filename = 'test_file_function' fid = open(filename + '.txt', 'w') start = time.mktime(time.strptime('2000', '%Y')) dt = 60 #One minute intervals t = 0.0 while t <= finaltime: t_string = time.strftime(time_format, time.gmtime(t+start)) fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600))) t += dt fid.close() #Convert ASCII file to NetCDF (Which is what we really like!) timefile2netcdf(filename) a = [0.0, 0.0] b = [4.0, 0.0] c = [0.0, 3.0] points = [a, b, c] vertices = [[0,1,2]] domain = Domain(points, vertices) #Check that domain.starttime isn't updated if later than file starttime but earlier #than file end time delta = 23 domain.starttime = start + delta F = file_function(filename + '.tms', domain, quantities = ['Attribute0', 'Attribute1', 'Attribute2']) assert allclose(domain.starttime, start+delta) #Now try interpolation with delta offset for i in range(20): t = i*10 q = F(t-delta) #Exact linear intpolation assert allclose(q[0], 2*t) if i%6 == 0: assert allclose(q[1], t**2) assert allclose(q[2], sin(t*pi/600)) #Check non-exact t = 90 #Halfway between 60 and 120 q = F(t-delta) assert allclose( (120**2 + 60**2)/2, q[1] ) assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] ) t = 100 #Two thirds of the way between between 60 and 120 q = F(t-delta) assert allclose( 2*120**2/3 + 60**2/3, q[1] ) assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] ) os.remove(filename + '.tms') os.remove(filename + '.txt') def test_xya_ascii(self): import time, os FN = 'xyatest' + str(time.time()) + '.xya' fid = open(FN, 'w') fid.write(' %s %s %s\n' %('a1', 'a2', 'a3') ) fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) ) fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) ) fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) ) fid.close() points, attributes = read_xya(FN, format = 'asc') assert allclose(points, [ [0,1], [1,0], [1,1] ]) assert allclose(attributes['a1'], [10,30,40.2]) assert allclose(attributes['a2'], [20,20,40.3]) assert allclose(attributes['a3'], [30,10,40.4]) os.remove(FN) def test_xya_ascii_w_names(self): import time, os FN = 'xyatest' + str(time.time()) + '.xya' fid = open(FN, 'w') fid.write(' %s %s %s\n' %('a1', 'a2', 'a3') ) fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) ) fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) ) fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) ) fid.close() points, attributes = read_xya(FN, format = 'asc') assert allclose(points, [ [0,1], [1,0], [1,1] ]) assert allclose(attributes['a1'], [10,30,40.2]) assert allclose(attributes['a2'], [20,20,40.3]) assert allclose(attributes['a3'], [30,10,40.4]) os.remove(FN) #Polygon stuff def test_polygon_function_constants(self): p1 = [[0,0], [10,0], [10,10], [0,10]] p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] f = Polygon_function( [(p1, 1.0)] ) z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1 assert allclose(z, [1,1,0,0]) f = Polygon_function( [(p2, 2.0)] ) z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2 assert allclose(z, [2,0,0,2]) #Combined f = Polygon_function( [(p1, 1.0), (p2, 2.0)] ) z = f([5, 5, 27, 35], [5, 9, 8, -5]) assert allclose(z, [2,1,0,2]) def test_polygon_function_callable(self): """Check that values passed into Polygon_function can be callable themselves. """ p1 = [[0,0], [10,0], [10,10], [0,10]] p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] f = Polygon_function( [(p1, test_function)] ) z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1 assert allclose(z, [10,14,0,0]) #Combined f = Polygon_function( [(p1, test_function), (p2, 2.0)] ) z = f([5, 5, 27, 35], [5, 9, 8, -5]) assert allclose(z, [2,14,0,2]) #Combined w default f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14) z = f([5, 5, 27, 35], [5, 9, 8, -5]) assert allclose(z, [2,14,3.14,2]) #Combined w default func f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = test_function) z = f([5, 5, 27, 35], [5, 9, 8, -5]) assert allclose(z, [2,14,35,2]) def test_point_on_line(self): #Endpoints first assert point_on_line( 0, 0, 0,0, 1,0 ) assert point_on_line( 1, 0, 0,0, 1,0 ) #Then points on line assert point_on_line( 0.5, 0, 0,0, 1,0 ) assert point_on_line( 0, 0.5, 0,1, 0,0 ) assert point_on_line( 1, 0.5, 1,1, 1,0 ) assert point_on_line( 0.5, 0.5, 0,0, 1,1 ) #Then points not on line assert not point_on_line( 0.5, 0, 0,1, 1,1 ) assert not point_on_line( 0, 0.5, 0,0, 1,1 ) #From real example that failed assert not point_on_line( 40,50, 40,20, 40,40 ) #From real example that failed assert not point_on_line( 40,19, 40,20, 40,40 ) def test_inside_polygon_main(self): #Simplest case: Polygon is the unit square polygon = [[0,0], [1,0], [1,1], [0,1]] assert inside_polygon( (0.5, 0.5), polygon ) assert not inside_polygon( (0.5, 1.5), polygon ) assert not inside_polygon( (0.5, -0.5), polygon ) assert not inside_polygon( (-0.5, 0.5), polygon ) assert not inside_polygon( (1.5, 0.5), polygon ) #Try point on borders assert inside_polygon( (1., 0.5), polygon, closed=True) assert inside_polygon( (0.5, 1), polygon, closed=True) assert inside_polygon( (0., 0.5), polygon, closed=True) assert inside_polygon( (0.5, 0.), polygon, closed=True) assert not inside_polygon( (0.5, 1), polygon, closed=False) assert not inside_polygon( (0., 0.5), polygon, closed=False) assert not inside_polygon( (0.5, 0.), polygon, closed=False) assert not inside_polygon( (1., 0.5), polygon, closed=False) #From real example (that failed) polygon = [[20,20], [40,20], [40,40], [20,40]] points = [ [40, 50] ] res = inside_polygon(points, polygon) assert len(res) == 0 polygon = [[20,20], [40,20], [40,40], [20,40]] points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ] res = inside_polygon(points, polygon) assert len(res) == 2 assert allclose(res, [0,1]) #More convoluted and non convex polygon polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] assert inside_polygon( (0.5, 0.5), polygon ) assert inside_polygon( (1, -0.5), polygon ) assert inside_polygon( (1.5, 0), polygon ) assert not inside_polygon( (0.5, 1.5), polygon ) assert not inside_polygon( (0.5, -0.5), polygon ) #Very convoluted polygon polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] assert inside_polygon( (5, 5), polygon ) assert inside_polygon( (17, 7), polygon ) assert inside_polygon( (27, 2), polygon ) assert inside_polygon( (35, -5), polygon ) assert not inside_polygon( (15, 7), polygon ) assert not inside_polygon( (24, 3), polygon ) assert not inside_polygon( (25, -10), polygon ) #Another combination (that failed) polygon = [[0,0], [10,0], [10,10], [0,10]] assert inside_polygon( (5, 5), polygon ) assert inside_polygon( (7, 7), polygon ) assert not inside_polygon( (-17, 7), polygon ) assert not inside_polygon( (7, 17), polygon ) assert not inside_polygon( (17, 7), polygon ) assert not inside_polygon( (27, 8), polygon ) assert not inside_polygon( (35, -5), polygon ) #Multiple polygons polygon = [[0,0], [1,0], [1,1], [0,1], [0,0], [10,10], [11,10], [11,11], [10,11], [10,10]] assert inside_polygon( (0.5, 0.5), polygon ) assert inside_polygon( (10.5, 10.5), polygon ) #FIXME: Fails if point is 5.5, 5.5 assert not inside_polygon( (0, 5.5), polygon ) #Polygon with a hole polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1], [0,0], [1,0], [1,1], [0,1], [0,0]] assert inside_polygon( (0, -0.5), polygon ) assert not inside_polygon( (0.5, 0.5), polygon ) def test_inside_polygon_vector_version(self): #Now try the vector formulation returning indices polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] res = inside_polygon( points, polygon, verbose=False ) assert allclose( res, [0,1,2] ) def test_outside_polygon(self): U = [[0,0], [1,0], [1,1], [0,1]] #Unit square assert not outside_polygon( [0.5, 0.5], U ) #evaluate to False as the point 0.5, 0.5 is inside the unit square assert outside_polygon( [1.5, 0.5], U ) #evaluate to True as the point 1.5, 0.5 is outside the unit square indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U ) assert allclose( indices, [1] ) #One more test of vector formulation returning indices polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] res = outside_polygon( points, polygon ) assert allclose( res, [3, 4] ) polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] res = outside_polygon( points, polygon ) assert allclose( res, [0, 4, 5] ) def test_outside_polygon2(self): U = [[0,0], [1,0], [1,1], [0,1]] #Unit square assert not outside_polygon( [0.5, 1.0], U, closed = True ) #evaluate to False as the point 0.5, 1.0 is inside the unit square assert outside_polygon( [0.5, 1.0], U, closed = False ) #evaluate to True as the point 0.5, 1.0 is outside the unit square def test_separate_points_by_polygon(self): U = [[0,0], [1,0], [1,1], [0,1]] #Unit square indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U ) assert allclose( indices, [0,2,1] ) assert count == 2 #One more test of vector formulation returning indices polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] res, count = separate_points_by_polygon( points, polygon ) assert allclose( res, [0,1,2,4,3] ) assert count == 3 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] res, count = separate_points_by_polygon( points, polygon ) assert allclose( res, [1,2,3,5,4,0] ) assert count == 3 def test_populate_polygon(self): polygon = [[0,0], [1,0], [1,1], [0,1]] points = populate_polygon(polygon, 5) assert len(points) == 5 for point in points: assert inside_polygon(point, polygon) #Very convoluted polygon polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] points = populate_polygon(polygon, 5) assert len(points) == 5 for point in points: assert inside_polygon(point, polygon) def test_populate_polygon_with_exclude(self): polygon = [[0,0], [1,0], [1,1], [0,1]] ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter points = populate_polygon(polygon, 5, exclude = [ex_poly]) assert len(points) == 5 for point in points: assert inside_polygon(point, polygon) assert not inside_polygon(point, ex_poly) #overlap polygon = [[0,0], [1,0], [1,1], [0,1]] ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]] points = populate_polygon(polygon, 5, exclude = [ex_poly]) assert len(points) == 5 for point in points: assert inside_polygon(point, polygon) assert not inside_polygon(point, ex_poly) #Multiple polygon = [[0,0], [1,0], [1,1], [0,1]] ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2]) assert len(points) == 20 for point in points: assert inside_polygon(point, polygon) assert not inside_polygon(point, ex_poly1) assert not inside_polygon(point, ex_poly2) #Very convoluted polygon polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]] points = populate_polygon(polygon, 20, exclude = [ex_poly]) assert len(points) == 20 for point in points: assert inside_polygon(point, polygon) assert not inside_polygon(point, ex_poly), '%s' %str(point) #------------------------------------------------------------- if __name__ == "__main__": suite = unittest.makeSuite(Test_Util,'test') #suite = unittest.makeSuite(Test_Util,'test_file_function_time') runner = unittest.TextTestRunner() runner.run(suite)