Changeset 4816
- Timestamp:
- Nov 15, 2007, 10:24:02 AM (17 years ago)
- Location:
- anuga_core/source/anuga/geospatial_data
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/geospatial_data/geospatial_data.py
r4776 r4816 11 11 size, shape 12 12 #from Array import tolist 13 from RandomArray import randint, seed 13 from RandomArray import randint, seed, get_seed 14 14 from copy import deepcopy 15 15 … … 670 670 671 671 # Set seed if provided, mainly important for unit test! 672 # plus recalcule seed when no seed provided. 672 673 if seed_num != None: 673 674 seed(seed_num,seed_num) 675 else: 676 seed() 677 if verbose: print "seed:", get_seed() 674 678 675 679 random_num = randint(0,self_size-1,(int(new_size),)) … … 1393 1397 alpha_list: the alpha values to test in a single list 1394 1398 1395 mesh_file: name of the created mesh file 1399 mesh_file: name of the created mesh file or if passed in will read it. 1400 NOTE, if there is a mesh file mesh_resolution, 1401 north_boundary, south... etc will be ignored. 1396 1402 1397 1403 mesh_resolution: the maximum area size for a triangle … … 1404 1410 1405 1411 USAGE: 1406 find_optimal_smoothing_parameter(data_file=fileName,1412 value, alpha = find_optimal_smoothing_parameter(data_file=fileName, 1407 1413 alpha_list=[0.0001, 0.01, 1], 1408 1414 mesh_file=None, … … 1430 1436 attribute_smoothed='elevation' 1431 1437 1432 if north_boundary is None or south_boundary is None or \1433 east_boundary is None or west_boundary is None:1434 no_boundary=True1435 else:1436 no_boundary=False1437 1438 if no_boundary is True:1439 msg= 'All boundaries must be defined'1440 raise msg1441 1442 1438 if mesh_file is None: 1443 1439 mesh_file='temp.msh' 1444 1440 1445 1446 1447 poly_topo = [[east_boundary,south_boundary], 1448 [east_boundary,north_boundary], 1449 [west_boundary,north_boundary], 1450 [west_boundary,south_boundary]] 1451 1452 create_mesh_from_regions(poly_topo, 1453 boundary_tags={'back': [2], 1454 'side': [1,3], 1455 'ocean': [0]}, 1456 maximum_triangle_area=mesh_resolution, 1457 filename=mesh_file, 1458 use_cache=cache, 1459 verbose=verbose) 1441 if north_boundary is None or south_boundary is None or \ 1442 east_boundary is None or west_boundary is None: 1443 no_boundary=True 1444 else: 1445 no_boundary=False 1446 1447 if no_boundary is True: 1448 msg= 'All boundaries must be defined' 1449 raise msg 1450 1451 poly_topo = [[east_boundary,south_boundary], 1452 [east_boundary,north_boundary], 1453 [west_boundary,north_boundary], 1454 [west_boundary,south_boundary]] 1455 1456 create_mesh_from_regions(poly_topo, 1457 boundary_tags={'back': [2], 1458 'side': [1,3], 1459 'ocean': [0]}, 1460 maximum_triangle_area=mesh_resolution, 1461 filename=mesh_file, 1462 use_cache=cache, 1463 verbose=verbose) 1464 1465 else: # if mesh file provided 1466 #test mesh file exists? 1467 if access(mesh_file,F_OK) == 0: 1468 msg="file %s doesn't exist!" %mesh_file 1469 raise IOError, msg 1460 1470 1461 1471 #split topo data … … 1530 1540 loglog(normal_cov_new[:,0],normal_cov_new[:,1]) 1531 1541 savefig(plot_name,dpi=300) 1532 1533 remove(mesh_file)1542 if mesh_file == 'temp.msh': 1543 remove(mesh_file) 1534 1544 1535 1545 return min(normal_cov_new[:,1]) , normal_cov_new[(argmin(normal_cov_new,axis=0))[1],0] -
anuga_core/source/anuga/geospatial_data/test_geospatial_data.py
r4745 r4816 1821 1821 assert allclose(A,[24, 18, 17, 11, 8]) 1822 1822 1823 1824 1825 # def test_find_optimal_smoothing_parameter1(self):1826 # from anuga.pmesh.mesh_interface import create_mesh_from_regions1827 # from anuga.shallow_water import Domain1828 # from anuga.shallow_water import Reflective_boundary1829 # from cmath import cos1830 #1831 # "make sphere domain"1832 #1833 # fileName = tempfile.mktemp(".csv")1834 # file = open(fileName,"w")1835 # file.write("x,y,elevation \n")1836 #1837 # for i in range(-5,6):1838 # for j in range(-5,6):1839 # z = abs(cos(((i*i) + (j*j))*.1)*2)1840 ## print 'x,y,f',i,j,z1841 # file.write("%s, %s, %s\n" %(i, j, z))1842 #1843 # file.close()1844 #1845 #1846 #1847 # poly = [[5,-5], [5,5], [-5,5], [-5,-5]]1848 # mesh_file='temp.msh'1849 #1850 # create_mesh_from_regions(poly,1851 # boundary_tags={'back': [2],1852 # 'side': [1,3],1853 # 'ocean': [0]},1854 # maximum_triangle_area=3,1855 # filename=mesh_file,1856 # use_cache=False,1857 # verbose=False)1858 #1859 # domain = Domain(mesh_file, use_cache=False, verbose=True)1860 # domain.set_name('test')1861 # domain.set_quantity('elevation', filename = fileName ,verbose=True, alpha=0.1)1862 #1863 #1864 # Br = Reflective_boundary(domain)1865 #1866 # domain.set_boundary({'back': Br,1867 # 'side': Br,1868 # 'ocean': Br})1869 #1870 # for t in domain.evolve(yieldstep = 1, finaltime = 1):1871 # domain.write_time()1872 1873 1823 def test_find_optimal_smoothing_parameter(self): 1874 1824 """ … … 1880 1830 from cmath import cos 1881 1831 1882 file Name = tempfile.mktemp(".csv")1883 file = open(file Name,"w")1832 filename = tempfile.mktemp(".csv") 1833 file = open(filename,"w") 1884 1834 file.write("x,y,elevation \n") 1885 1835 … … 1893 1843 file.close() 1894 1844 1895 value, alpha = find_optimal_smoothing_parameter(data_file=file Name,1845 value, alpha = find_optimal_smoothing_parameter(data_file=filename, 1896 1846 alpha_list=[0.0001, 0.01, 1], 1897 1847 mesh_file=None, … … 1905 1855 verbose=False) 1906 1856 1907 1908 results = Geospatial_data(fileName) 1909 os.remove(fileName) 1910 1911 # print value, alpha 1857 os.remove(filename) 1858 1859 # print value, alpha 1912 1860 assert (alpha==0.01) 1913 1861 1862 def test_find_optimal_smoothing_parameter1(self): 1863 """ 1864 Creates a elevation file represting hill (sort of) and 1865 Then creates a mesh file and passes the mesh file and the elevation 1866 file to find_optimal_smoothing_parameter for 3 different alphas, 1867 1868 NOTE the random number seed is provided to control the results 1869 """ 1870 from cmath import cos 1871 from anuga.pmesh.mesh_interface import create_mesh_from_regions 1872 1873 filename = tempfile.mktemp(".csv") 1874 file = open(filename,"w") 1875 file.write("x,y,elevation \n") 1876 1877 for i in range(-5,6): 1878 for j in range(-5,6): 1879 #this equation made surface like a circle ripple 1880 z = abs(cos(((i*i) + (j*j))*.1)*2) 1881 # print 'x,y,f',i,j,z 1882 file.write("%s, %s, %s\n" %(i, j, z)) 1883 1884 file.close() 1885 poly=[[5,5],[5,-5],[-5,-5],[-5,5]] 1886 internal_poly=[[[[1,1],[1,-1],[-1,-1],[-1,1]],.5]] 1887 mesh_filename= tempfile.mktemp(".msh") 1888 1889 create_mesh_from_regions(poly, 1890 boundary_tags={'back': [2], 1891 'side': [1,3], 1892 'ocean': [0]}, 1893 maximum_triangle_area=3, 1894 interior_regions=internal_poly, 1895 filename=mesh_filename, 1896 use_cache=False, 1897 verbose=False) 1898 1899 value, alpha = find_optimal_smoothing_parameter(data_file=filename, 1900 alpha_list=[0.0001, 0.01, 1], 1901 mesh_file=mesh_filename, 1902 plot_name=None, 1903 seed_num=174, 1904 verbose=False) 1905 1906 os.remove(filename) 1907 os.remove(mesh_filename) 1908 1909 # print value, alpha 1910 assert (alpha==0.01) 1911 1912 def test_find_optimal_smoothing_parameter2(self): 1913 """ 1914 Tests requirement that mesh file must exist or IOError is thrown 1915 1916 NOTE the random number seed is provided to control the results 1917 """ 1918 from cmath import cos 1919 from anuga.pmesh.mesh_interface import create_mesh_from_regions 1920 1921 filename = tempfile.mktemp(".csv") 1922 mesh_filename= tempfile.mktemp(".msh") 1923 1924 try: 1925 value, alpha = find_optimal_smoothing_parameter(data_file=filename, 1926 alpha_list=[0.0001, 0.01, 1], 1927 mesh_file=mesh_filename, 1928 plot_name=None, 1929 seed_num=174, 1930 verbose=False) 1931 except IOError: 1932 pass 1933 else: 1934 self.failUnless(0 ==1, 'Error not thrown error!') 1914 1935 1915 1936 … … 1917 1938 1918 1939 #suite = unittest.makeSuite(Test_Geospatial_data, 'test_write_csv_attributes_lat_long') 1919 suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')1940 # suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter') 1920 1941 # suite = unittest.makeSuite(Test_Geospatial_data, 'test_split') 1921 #suite = unittest.makeSuite(Test_Geospatial_data, 'test')1942 suite = unittest.makeSuite(Test_Geospatial_data, 'test') 1922 1943 runner = unittest.TextTestRunner() #verbosity=2) 1923 1944 runner.run(suite)
Note: See TracChangeset
for help on using the changeset viewer.