Changeset 4091

Timestamp:

Dec 19, 2006, 9:19:27 AM (18 years ago)

Author:

nick

Message:

updates to dampier and perth

Location:

anuga_work/production

Files:

• dampier_2006/project.py (modified) (2 diffs)
• perth_2006/build_perth.py (modified) (7 diffs)
• perth_2006/project.py (modified) (9 diffs)
• perth_2006/run_perth.py (modified) (8 diffs)

anuga_work/production/dampier_2006/project.py

@@ -203 +203 @@
 n_max_area = 7725000
 
-e_min_area_mom = 442200
-e_max_area_mom = 497028
-n_min_area_mom = 7720095
-n_max_area_mom = 7750500
-
 poly_facility = read_polygon(polygons_dir+'facility.csv')
 
@@ -224 +219 @@
 clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs')
 
-plot_polygons([bounding_polygon,poly_facility,poly_pipeline,poly_interior,poly_coast],'polys',verbose=True)
 #Interior regions
 karratha_south = degminsec2decimal_degrees(-20,44,0)

anuga_work/production/perth_2006/build_perth.py

@@ -57 +57 @@
 print"project.bounding_polygon",project.bounding_polygon
 print"project.combined_dir_name",project.combined_dir_name
-'''
+
 # topography directory filenames
 onshore_in_dir_name = project.onshore_in_dir_name
@@ -66 +66 @@
 island_in_dir_name3 = project.island_in_dir_name3
 offshore_in_dir_name = project.offshore_in_dir_name
+offshore1_in_dir_name = project.offshore1_in_dir_name
 
 onshore_dir_name = project.onshore_dir_name
@@ -76 +77 @@
 
 # creates DEM from asc data
+print "creates DEMs from asc data"
 convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
 convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
@@ -83 +85 @@
 
 #creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
 dem2pts(onshore_dir_name,
 #        easting_min=project.eastingmin,
@@ -97 +100 @@
 dem2pts(island_dir_name3, use_cache=True, verbose=True)
 
-print'create Geospatial data objects from topographies'
+print'create Geospatial data1 objects from topographies'
 G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
+print'create Geospatial data2 objects from topographies'
 G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
+print'create Geospatial data3 objects from topographies'
 G3 = Geospatial_data(file_name = island_dir_name + '.pts')
+print'create Geospatial data4 objects from topographies'
 G4 = Geospatial_data(file_name = island_dir_name1 + '.pts')
+print'create Geospatial data5 objects from topographies'
 G5 = Geospatial_data(file_name = island_dir_name2 + '.pts')
+print'create Geospatial data6 objects from topographies'
 G6 = Geospatial_data(file_name = island_dir_name3 + '.pts')
+print'create Geospatial data7 objects from topographies'
 G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
+print'create Geospatial data8 objects from topographies'
+G_off1 = Geospatial_data(file_name = offshore1_in_dir_name + '.xya')
 
 print'add all geospatial objects'
-G = G1 + G2 + G3 + G4 + G5 + G6 + G_off 
+G = G1 + G2 + G3 + G4 + G5 + G6 + G_off + G_off1 
 
 print'clip combined geospatial object by bounding polygon'
@@ -117 +128 @@
 if access(project.topographies_dir,F_OK) == 0:
     mkdir (project.topographies_dir)
-#G_clipped.export_points_file(project.combined_dir_name + '.pts')
-G_clipped.export_points_file(project.combined_dir_name + '.xya')
+G_clipped.export_points_file(project.combined_dir_name + '.pts')
+#G_clipped.export_points_file(project.combined_dir_name + '.xya')
 
 '''
@@ -130 +141 @@
 G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
 
-'''
+
 #-------------------------------------------------------------------------
 # Convert URS to SWW file for boundary conditions 

anuga_work/production/perth_2006/project.py

@@ -46 +46 @@
 coast_name = 'waterline'
 offshore_name = 'perth_bathymetry'
+offshore1_name = 'missing_fairsheets'
 
 #final topo name
 combined_name ='perth_combined_elevation'
+combined_smaller_name = 'perth_combined_elevation_smaller'
+
 
 topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
@@ -63 +66 @@
 coast_in_dir_name = topographies_in_dir + coast_name
 offshore_in_dir_name = topographies_in_dir + offshore_name
+offshore1_in_dir_name = topographies_in_dir + offshore1_name
 
 onshore_dir_name = topographies_dir + onshore_name
@@ -76 +80 @@
 combined_dir_name = topographies_dir + combined_name
 combined_time_dir_name = topographies_time_dir + combined_name
+combined_smaller_name_dir = topographies_dir + combined_smaller_name
 #combined_time_dir_final_name = topographies_time_dir + combined_final_name
 
 
-#Derive subdirectories and filenames
 
 meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep
@@ -88 +92 @@
 
 
-boundaries_source = 'mag_9_corrected'
+boundaries_source = '????'
 #boundaries locations
 boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
@@ -96 +100 @@
 #boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
 #boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
-#ideas
-#boundaries_time_dir = boundaries_in_dir+'urs'+sep+boundaries_source+sep
 
 #output locations
@@ -111 +113 @@
 
 
-
-
-
-
-
-
 ###############################
 # Domain definitions
@@ -123 +119 @@
 
 bounding_polygon = read_polygon(polygons_dir+'bounding_poly.csv')
+res_bounding = 1000000
 
 refzone = 50 
 
 
-#Interior regions
+###############################
+# Interior region definitions
+###############################
 
-#cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
+poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20pts.csv')
+res_pos20_neg20 = 500000
+
+poly_cbd = read_polygon(polygons_dir+'cbd_pts.csv')
+res_cbd = 50000
+
+poly_penguin = read_polygon(polygons_dir+'penguin_pts.csv')
+res_penguin = 50000
 #assert zone == refzone
 
-poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20pts.csv')
+trigs_bound = polygon_area(bounding_polygon)/res_bounding
+trigs_pos = polygon_area(poly_pos20_neg20)/res_pos20_neg20
+trigs_cbd = polygon_area(poly_cbd)/res_cbd
+trigs_penguin = polygon_area(poly_penguin)/res_penguin
+trigs_min = trigs_bound + trigs_pos + trigs_cbd + trigs_penguin
 
-poly_cbd = read_polygon(polygons_dir+'cbd_pts.csv')
-
-poly_penguin = read_polygon(polygons_dir+'penguin_pts.csv')
-#assert zone == refzone
-
-print "bounding_polygon", bounding_polygon
-print 'Area of bounding poly', polygon_area(bounding_polygon)/1000000.0
-print 'Area of pos20_neg20pts', polygon_area(poly_pos20_neg20)/1000000.0
-print 'Area of poly_cbd', polygon_area(poly_cbd)/1000000.0
-print 'Area of poly_penguin', polygon_area(poly_penguin)/1000000.0
+print 'Area of bounding poly', trigs_bound
+print 'Area of pos20_neg20pts', trigs_pos
+print 'Area of poly_cbd', trigs_cbd
+print 'Area of poly_penguin', trigs_penguin
+print 'min number triangles', trigs_min
 
 
@@ -149 +154 @@
 # Clipping regions for export to asc and regions for clipping data
 ###################################################################
-'''
+
 # exporting asc grid
 eastingmin = 406215.87
@@ -156 +161 @@
 northingmax = 8032834.52
 
-###############################
-# Interior region definitions
-###############################
 
-# perth digitized polygons
-poly_perth1 = read_polygon(polygondir+'perth_Local_Polygon_update.csv')
-poly_perth2 = read_polygon(polygondir+'perth_Close2_update.csv')
-poly_perth3 = read_polygon(polygondir+'perth_Coast_update.csv')
-#poly_perth4 = read_polygon(polygondir+'Cable_Beach_revised.csv')
 
-plot_polygons([polyAll,poly_perth1,poly_perth2,poly_perth3],'boundingpoly2',verbose=False)
-print 'Area of local polygon', polygon_area(poly_perth1)/1000000.0
-print 'Area of close polygon', polygon_area(poly_perth2)/1000000.0
-print 'Area of coastal polygon', polygon_area(poly_perth3)/1000000.0
-#print 'Area of cable beach polygon', polygon_area(poly_perth4)/1000000.0
-'''
-

anuga_work/production/perth_2006/run_perth.py

@@ -31 +31 @@
 
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
-
 from anuga.geospatial_data.geospatial_data import *
 from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
 from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
+from anuga.utilities.polygon import plot_polygons, polygon_area
 
 # Application specific imports
@@ -47 +47 @@
 # filenames
 
-build_time = '20061107_063840_build'
-bound_time = '20061102_221245_build'
-
 boundaries_name = project.scenario
 meshes_dir_name = project.meshes_dir_name+'.msh'
-#source_dir = project.boundarydir
 #boundaries_time_dir_name = project.boundaries_time_dir_name
 boundaries_dir_name = project.boundaries_dir_name
@@ -67 +63 @@
 
 print 'USER: ', project.user
+print 'min triangles', project.trigs_min,
+print 'Note: This is generally about 20% less than the final amount'
 
 #--------------------------------------------------------------------------
@@ -76 +74 @@
 
 if myid == 0:
+
     print 'start create mesh from regions'
-    interior_regions = [[project.poly_pos20_neg20, 500000],
-                    [project.poly_cbd, 50000],
-                    [project.poly_penguin, 50000]]    
+    interior_regions = [[project.poly_pos20_neg20, project.res_pos20_neg20],
+                    [project.poly_cbd, project.res_cbd],
+                    [project.poly_penguin, project.res_penguin]]    
 #                    [project.poly_interior, 1000]]       
 
@@ -89 +88 @@
 #               verbose = True)
 #    import sys; sys.exit()
+
     print 'start create mesh from regions'
-
     create_mesh_from_regions(project.bounding_polygon,
                          boundary_tags={'back': [4], 'side': [0,3],
                                         'ocean': [1, 2]},
-                         maximum_triangle_area=1000000,
+                         maximum_triangle_area=project.res_bounding,
                          interior_regions=interior_regions,
 #                         filename=meshes_time_dir_name,
@@ -157 +156 @@
 
     domain.set_quantity('elevation', 
-                    filename = project.combined_dir_name+'.xya',
+#                    filename = project.combined_dir_name+'.xya',
+                    filename = project.combined_smaller_name_dir+'.xya',
                     use_cache = True,
                     verbose = True,
@@ -201 +201 @@
 domain.set_boundary({'back': Br,
                      'side': Bd,
-                     'ocean': Bf}) 
+                     'ocean': Bd}) 
 print'finish set boundary'
 
@@ -214 +214 @@
     domain.write_boundary_statistics(tags = 'ocean')
     if allclose(t, 120):
-        domain.set_boundary({'back': Br, 'side': Bo, 'ocean': Bo})
+        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
 
     if allclose(t, 1020):
-        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
+        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})