Index: /anuga_work/production/hobart_2006/project.py
===================================================================
--- /anuga_work/production/hobart_2006/project.py	(revision 3700)
+++ /anuga_work/production/hobart_2006/project.py	(revision 3701)
@@ -26,36 +26,40 @@
 scenario_dir_name = 'hobart_tsunami_scenario_2006'
 
-# data provided by Tas SES and checked by NM&I
-#onshore_name = 'hob3_topo' # original
-onshore_name = 'hob5_topo' # 12.5m grid and clipped to 100m elevation or 3000m from coast
-#onshore_name_25 = 'hob5_topo_25m' # 25m grid and clipped to 100m elevation or 3000m from coast
-onshore_name_25 = 'hob6_topo_25m' # 25m grid NOT clipped
-offshore_name_tas1 = 'derwent_2m'
-offshore_name_tas2 = 'derwent_5m'
-offshore_name_tas3 = 'south_east_tas' #actually this is AHO
-offshore_name_tas4 = 'hobart_1m'
-
-# AHO data and checked by NM&I
-offshore_name1 = 'xy100003760'
-offshore_name2 = 'xy100003761'
-offshore_name3 = 'xy100003762'
-offshore_name4 = 'xy100003907'
-offshore_name5 = 'xy100003908'
-offshore_name6 = 'xy100003909'
-offshore_name7 = 'xy100003910'
-offshore_name8 = 'xy100003932'
-offshore_name9 = 'xy100003933'
-offshore_name10 = 'xy100003934'
-offshore_name11 = 'xy100003935'
-offshore_name12 = 'xy100003936'
-offshore_name13 = 'xy100003964'
-offshore_name14 = 'xy100014250'
-offshore_name15 = 'xy100014253'
-offshore_name16 = 'xy100016142'
-
-# developed by NM&I
-coast_name = 'coastline_points'
-
-boundary_basename = 'puysegur'
+
+##
+### data provided by Tas SES and checked by NM&I
+###onshore_name = 'hob3_topo' # original
+###onshore_name_25 = 'hob5_topo_25m' # 25m grid and clipped to 100m elevation or 3000m from coast
+##onshore_name = 'hob5_topo' # 12.5m grid and clipped to 100m elevation or 3000m from coast
+##onshore_name_25 = 'hob6_topo_25m' # 25m grid NOT clipped
+##offshore_name_tas1 = 'derwent_2m'
+##offshore_name_tas2 = 'derwent_5m'
+##offshore_name_tas3 = 'south_east_tas' #actually this is AHO
+##offshore_name_tas4 = 'hobart_1m'
+##
+### AHO data and checked by NM&I
+##offshore_name1 = 'xy100003760'
+##offshore_name2 = 'xy100003761'
+##offshore_name3 = 'xy100003762'
+##offshore_name4 = 'xy100003907'
+##offshore_name5 = 'xy100003908'
+##offshore_name6 = 'xy100003909'
+##offshore_name7 = 'xy100003910'
+##offshore_name8 = 'xy100003932'
+##offshore_name9 = 'xy100003933'
+##offshore_name10 = 'xy100003934'
+##offshore_name11 = 'xy100003935'
+##offshore_name12 = 'xy100003936'
+##offshore_name13 = 'xy100003964'
+##offshore_name14 = 'xy100014250'
+##offshore_name15 = 'xy100014253'
+##offshore_name16 = 'xy100016142'
+##
+### developed by NM&I
+##coast_name = 'coastline_points'
+
+boundary_basename = 'puysegur' # Mw 8.7
+#boundary_basename = 'puysegur_clip' # Mw 8.5
+
 
 #swollen/ all data output
@@ -76,8 +80,8 @@
 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
 
-gauge_filename = gaugedir + 'hobart_gauges.csv'
-buildings_filename = gaugedir + 'hobart_res.csv'
-buildings_filename_damage_out = 'hobart_res_modified.csv'
-gaugetimeseries = gaugedir + 'hobart'
+gauge_filename = gaugedir + 'hobart_gauges_final.csv'
+#buildings_filename = gaugedir + 'hobart_res.csv'
+#buildings_filename_damage_out = 'hobart_res_modified.csv'
+#gaugetimeseries = gaugedir + 'hobart'
 
 # boundary source data
@@ -88,38 +92,43 @@
 meshname = outputtimedir + 'mesh_' + basename
 
-onshore_dem_name = datadir + onshore_name
-onshore_dem_name_25 = datadir + onshore_name_25
-all_onshore_dem_name = datadir + 'combined_onshore'
-offshore_dem_name_local1 = datadir + offshore_name_tas1
-offshore_dem_name_local2 = datadir + offshore_name_tas2
-offshore_dem_name_local3 = datadir + offshore_name_tas3
-offshore_dem_name_local4 = datadir + offshore_name_tas4
-offshore_dem_name_aho1 = datadir + offshore_name1
-offshore_dem_name_aho2 = datadir + offshore_name2
-offshore_dem_name_aho3 = datadir + offshore_name3
-offshore_dem_name_aho4 = datadir + offshore_name4
-offshore_dem_name_aho5 = datadir + offshore_name5
-offshore_dem_name_aho6 = datadir + offshore_name6
-offshore_dem_name_aho7 = datadir + offshore_name7
-offshore_dem_name_aho8 = datadir + offshore_name8
-offshore_dem_name_aho9 = datadir + offshore_name9
-offshore_dem_name_aho10 = datadir + offshore_name10
-offshore_dem_name_aho11 = datadir + offshore_name11
-offshore_dem_name_aho12 = datadir + offshore_name12
-offshore_dem_name_aho13 = datadir + offshore_name13
-offshore_dem_name_aho14 = datadir + offshore_name14
-offshore_dem_name_aho15 = datadir + offshore_name15
-offshore_dem_name_aho16 = datadir + offshore_name16
-coast_dem_name = datadir + coast_name
+# Necessary if using point datasets, rather than grid
+##onshore_dem_name = datadir + onshore_name
+##onshore_dem_name_25 = datadir + onshore_name_25
+##all_onshore_dem_name = datadir + 'combined_onshore'
+##offshore_dem_name_local1 = datadir + offshore_name_tas1
+##offshore_dem_name_local2 = datadir + offshore_name_tas2
+##offshore_dem_name_local3 = datadir + offshore_name_tas3
+##offshore_dem_name_local4 = datadir + offshore_name_tas4
+##offshore_dem_name_aho1 = datadir + offshore_name1
+##offshore_dem_name_aho2 = datadir + offshore_name2
+##offshore_dem_name_aho3 = datadir + offshore_name3
+##offshore_dem_name_aho4 = datadir + offshore_name4
+##offshore_dem_name_aho5 = datadir + offshore_name5
+##offshore_dem_name_aho6 = datadir + offshore_name6
+##offshore_dem_name_aho7 = datadir + offshore_name7
+##offshore_dem_name_aho8 = datadir + offshore_name8
+##offshore_dem_name_aho9 = datadir + offshore_name9
+##offshore_dem_name_aho10 = datadir + offshore_name10
+##offshore_dem_name_aho11 = datadir + offshore_name11
+##offshore_dem_name_aho12 = datadir + offshore_name12
+##offshore_dem_name_aho13 = datadir + offshore_name13
+##offshore_dem_name_aho14 = datadir + offshore_name14
+##offshore_dem_name_aho15 = datadir + offshore_name15
+##offshore_dem_name_aho16 = datadir + offshore_name16
+##coast_dem_name = datadir + coast_name
+
+
 # addition once total grid delivered
 onshore_offshore_dem_name_25 = datadir + '25m_se_tas' #25m grid
-onshore_offshore_dem_name = datadir + '50m_se_tas' #50m grid
+onshore_offshore_dem_name    = datadir + '50m_se_tas' #50m grid
+
 # output names
-bruny_dem_name_25 = datadir + 'bruny_25_dem'
-hobart_dem_name_25 = datadir + 'hobart_25_dem'
-combined_dem_name = datadir + 'hobart_combined_elevation'
+bruny_dem_name_25   = datadir + 'bruny_25_dem'
+hobart_dem_name_25  = datadir + 'hobart_25_dem'
+combined_dem_name   = datadir + 'hobart_combined_elevation'
 combined_dem_name_2 = datadir + 'hobart_combined_elevation_2'
 
-outputname = outputtimedir + basename  #Used by post processing
+#outputname = outputtimedir + basename  #Used by post processing
+
 
 ###############################
@@ -127,11 +136,11 @@
 ###############################
 
-# bounding box
+# bounding box for clipping MOST output (much bigger than study area)
 south = degminsec2decimal_degrees(-44,45,0)
 north = degminsec2decimal_degrees(-42,0,0)
-west = degminsec2decimal_degrees(146,45,0)
-east = degminsec2decimal_degrees(148,25,0) #degminsec2decimal_degrees(148,15,0)
-
-#Main Domain of Hobart: first run JS 18/9/06
+west  = degminsec2decimal_degrees(146,45,0)
+east  = degminsec2decimal_degrees(148,25,0) 
+
+###Main Domain of Hobart: 
 d0 = [south, west]
 d1 = [south, east]
@@ -139,5 +148,4 @@
 d3 = [north, west]
 polyAll2, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
-print "polyAll2", polyAll2
 refzone = zone
 
@@ -145,7 +153,4 @@
 #polyAll = [[520000, 5170000],[580000, 5170000],[580000, 5200000],[590000,5240000],[520000,5260000]]
 
-# Third run - afternoon Wed 27 Sep; surrounds -100m and 20mish elevation
-#polyAll = read_polygon(polygondir+'new_extent.csv')
-#plot_polygons([polyAll],'boundingpoly',verbose=False)
 
 # Mark run - morning Fri 29 Sep; surrounds -100m and 20mish elevation
@@ -158,28 +163,10 @@
 ###################################################################
 
-# region to export for inundation map
-e_min_area = 520000#490000
-e_max_area = 580000#580000
-n_min_area = 5180000#5160000
-n_max_area = 5260000#5275000
-
-# clipping 12.5m onshore data set
-eastingmin = 520000
-eastingmax = 536000
-northingmin = 5245000
-northingmax = 5265000
-
-### clipping 25m data set - Hobart
-##eastingmin25  = 522000
-##eastingmax25  = 555000
-##northingmin25 = 5234000
-##northingmax25 = 5260000
 
 # clipping 25m data set - Hobart
-eastingmin25  = 524208.387
-eastingmax25  = 554867.24
+eastingmin25 = 524208.387
+eastingmax25 = 554867.24
 northingmin25 = 5229154.555
 northingmax25 = 5258511.857
-
 
 # clipping 25m data set - Bruny
@@ -225,5 +212,4 @@
 poly_hobart2 = read_polygon(polygondir+'Hob_poly2.csv')
 poly_hobart3 = read_polygon(polygondir+'Hob_poly3.csv')
-poly_hobart4 = read_polygon(polygondir+'Hob_poly4.csv')
-
-plot_polygons([polyAll, poly_hobart1,poly_hobart2,poly_hobart3, poly_hobart4],'boundingpoly2',verbose=False)
+
+plot_polygons([polyAll, poly_hobart1,poly_hobart2,poly_hobart3],'boundingpoly2',verbose=False)
Index: /anuga_work/production/hobart_2006/run_hobart.py
===================================================================
--- /anuga_work/production/hobart_2006/run_hobart.py	(revision 3700)
+++ /anuga_work/production/hobart_2006/run_hobart.py	(revision 3701)
@@ -10,5 +10,6 @@
 Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
 """
-#-------------------------------------------------------------------------------# Import necessary modules
+#-------------------------------------------------------------------------------
+# Import necessary modules
 #-------------------------------------------------------------------------------
 
@@ -49,5 +50,4 @@
 #used to catch screen output to file
 sys.stdout = Screen_Catcher(screen_output_name)
-#sys.stderr = Screen_Catcher(screen_output_name)
 sys.stderr = Screen_Catcher(screen_error_name)
 
@@ -76,4 +76,5 @@
 # 25m data (clipping the around the Hobart area)
 convert_dem_from_ascii2netcdf(onshore_offshore_dem_name_25, use_cache=True, verbose=True)
+
 # creates pts file for 25m data around Hobart
 dem2pts(onshore_offshore_dem_name_25, project.hobart_dem_name_25,
@@ -91,5 +92,6 @@
 
 # 25m data (clipping the around site 24 on Bruny Island)
-convert_dem_from_ascii2netcdf(onshore_offshore_dem_name_25, use_cache=True, verbose=True)
+#convert_dem_from_ascii2netcdf(onshore_offshore_dem_name_25, use_cache=True, verbose=True)
+
 # creates pts file for 25m data around site 24 at Bruny Island
 dem2pts(onshore_offshore_dem_name_25, project.bruny_dem_name_25,
@@ -106,7 +108,9 @@
                                  project.combined_dem_name_2 + '.pts')
 
+# If working with points, then add separate datasets together here
 # create geospatial data set and export 
 #G = Geospatial_data(file_name = project.onshore_offshore_dem_name + '.pts')
 #G.export_points_file(project.combined_dem_name + '.pts')
+
 
 #----------------------------------------------------------------------------
@@ -119,14 +123,18 @@
 
 # use 75 for onshore components (12.5m DEM)
-island_res = 35000
-hobart_res = 5000
-peninsula_res = 35000
+# Resolution refers to max triangle area in a region
+#island_res    = 35000
+hobart_res    =  7000
+#peninsula_res = 35000
 interior_regions = [[project.poly_hobart1, hobart_res],
                     [project.poly_hobart2, hobart_res],
-                    [project.poly_hobart3, hobart_res],
-                    [project.poly_hobart4, hobart_res]]
+                    [project.poly_hobart3, hobart_res]]
 
 print 'number of interior regions', len(interior_regions)
 
+
+# Number tags correspond to number segments from bounding polygon
+# maximum_triangle_area refers to max triangle area in a region
+# Interion_regions either comment out (perhaps empty list)
 from caching import cache
 _ = cache(create_mesh_from_regions,
@@ -138,5 +146,5 @@
 			      'e12': [12], 'e13': [13], 'e14': [14],
 			      'e15': [15]},
-           'maximum_triangle_area': 200000,
+           'maximum_triangle_area': 750000,
            'filename': meshname,
            'interior_regions': interior_regions},
@@ -159,4 +167,5 @@
 domain.set_store_vertices_uniquely(False)  # for writting to sww
 
+
 #-------------------------------------------------------------------------------                                 
 # Setup initial conditions
@@ -164,12 +173,9 @@
 
 tide = 0.0
-
 domain.set_quantity('stage', tide)
 domain.set_quantity('friction', 0.0) 
 
 domain.set_quantity('elevation', 
-#                    filename = project.onshore_dem_name + '.pts',
                     filename = project.combined_dem_name_2 + '.pts',
-#                    filename = project.offshore_dem_name + '.pts',
                     use_cache = True,
                     verbose = True,
@@ -177,4 +183,5 @@
                     )
 
+
 #-------------------------------------------------------------------------------                                 
 # Setup boundary conditions
@@ -185,25 +192,25 @@
 south = project.south 
 north = project.north
-west = project.west
-east = project.east
-
-#note only need to do when an SWW file for the MOST boundary doesn't exist
-cache(ferret2sww,
-      (source_dir + project.boundary_basename,
-       source_dir + project.boundary_basename), 
-      {'verbose': True,
-       'minlat': south,
-       'maxlat': north,
-       'minlon': west,
-       'maxlon': east,
-#       'origin': project.mesh_origin,
-       'origin': domain.geo_reference.get_origin(),
-       'mean_stage': tide,
-       'zscale': 1,                 #Enhance tsunami
-       'fail_on_NaN': False,
-       'inverted_bathymetry': True},
-      #evaluate = True,
-       verbose = True,
-       dependencies = source_dir + project.boundary_basename + '.sww')
+west  = project.west
+east  = project.east
+
+###note only need to do when an SWW file for the MOST boundary doesn't exist
+##cache(ferret2sww,
+##      (source_dir + project.boundary_basename,
+##       source_dir + project.boundary_basename), 
+##      {'verbose': True,
+##       'minlat': south,
+##       'maxlat': north,
+##       'minlon': west,
+##       'maxlon': east,
+###       'origin': project.mesh_origin,
+##       'origin': domain.geo_reference.get_origin(),
+##       'mean_stage': tide,
+##       'zscale': 1,                 #Enhance tsunami
+##       'fail_on_NaN': False,
+##       'inverted_bathymetry': True},
+##      #evaluate = True,
+##       verbose = True,
+##       dependencies = source_dir + project.boundary_basename + '.sww')
 
 
@@ -220,15 +227,9 @@
                    f=lambda t: [(60<t<480)*10, 0, 0])
 
-# for MOST BC
-#domain.set_boundary( {'top': Bd, 'left': Bd,
-#                      'bottom': Bf, 'right': Bf} )
-
-# for testing
-#domain.set_boundary( {'topr': Bd, 'left': Bd, 'top': Bd,
-#                      'bottom': Bw, 'bright': Bd} )
-domain.set_boundary( {'e0': Bd, 'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
-			'e5': Bd, 'e6': Bd, 'e7': Bd, 'e8': Bd, 'e9': Bd,
-			'e10': Bd, 'e11': Bd, 'e12': Bf, 'e13': Bf, 'e14': Bf,
-			'e15': Bf} )
+domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
+		      'e5': Bd,  'e6': Bd, 'e7': Bd, 'e8': Bd, 'e9': Bd,
+		      'e10': Bd, 'e11': Bd, 'e12': Bf, 'e13': Bf, 'e14': Bf,
+ 		      'e15': Bf} )
+
 
 #-------------------------------------------------------------------------------                                 
@@ -240,15 +241,15 @@
 for t in domain.evolve(yieldstep = 240, finaltime = 6800): 
     domain.write_time()
-    domain.write_boundary_statistics(tags = 'bottom')     
+    domain.write_boundary_statistics(tags = 'e14')     
 
 for t in domain.evolve(yieldstep = 30, finaltime = 9000
                        ,skip_initial_step = True): 
     domain.write_time()
-    domain.write_boundary_statistics(tags = 'bottom')     
+    domain.write_boundary_statistics(tags = 'e14')     
 
 for t in domain.evolve(yieldstep = 240, finaltime = 15000
                        ,skip_initial_step = True): 
     domain.write_time()
-    domain.write_boundary_statistics(tags = 'bottom') 
+    domain.write_boundary_statistics(tags = 'e14') 
     
 print 'That took %.2f seconds' %(time.time()-t0)
