Index: /anuga_work/production/broome_2006/project.py
===================================================================
--- /anuga_work/production/broome_2006/project.py	(revision 3787)
+++ /anuga_work/production/broome_2006/project.py	(revision 3788)
@@ -1,13 +1,24 @@
+# -*- coding: cp1252 -*-
 """Common filenames and locations for topographic data, meshes and outputs.
 """
 
 from os import sep, environ, getenv, getcwd
-from os.path import expanduser, basename
-#from anuga.utilities.polygon import read_polygon
+from os.path import expanduser
 import sys
-from anuga.pmesh.create_mesh import convert_from_latlon_to_utm
-from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
-from time import localtime, strftime
-from anuga.geospatial_data.geospatial_data import *
+from time import localtime, strftime, gmtime
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
+
+if sys.platform == 'win32':
+    home = getenv('INUNDATIONHOME')
+    user = getenv('USERPROFILE')
+
+else:    
+    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
+    user = getenv('LOGNAME')
+    print 'USER:', user
+
+# INUNDATIONHOME is the inundation directory, not the data directory.
+home += sep +'data'
 
 #Making assumptions about the location of scenario data
@@ -15,19 +26,17 @@
 scenario_dir_name = 'broome_tsunami_scenario_2006'
 
-# all data to be delivered by National Mapping
-# onshore data from 30m DTED level 2
-onshore_name_dted = 'broome_onshore_30m_dted'  
-onshore_name_dli = 'broome_onshore_20m_dli' 
+# onshore data provided by WA DLI
+onshore_name = '' # original
 
-# offshore data from GA digitised charts
-offshore_name1 = 'broome_offshore_points'
+# offshore data provided by WA DPI
+offshore_name_dpi1 = ''
 
-# offshore data from AHO fairsheets
-offshore_name2 = 'broome_offshore_points_fairsheet'
+# AHO data
+offshore_name1 = 'xy100003760'
 
-# coastline developed from aerial photography and 1.5m DLI contour
-coast_name = 'broome_coastline_points'
+# developed by NM&I
+coast_name = 'coastline_points'
 
-boundary_basename = 'SU-AU_clip'
+boundary_basename = 'SU-AU' # Mw ?
 
 #swollen/ all data output
@@ -35,126 +44,81 @@
 codename = 'project.py'
 
-if sys.platform == 'win32':
-    home = getenv('INUNDATIONHOME') #Sandpit's parent dir     
-    user = getenv('USERPROFILE')
-else:
-    # update to perlite 2
-    home = getenv('INUNDATIONHOME', sep+'d'+sep+'cit'+sep+'2'+sep+'cit'+sep+'inundation'+sep+'data')     
-    user = getenv('LOGNAME')
-
-# INUNDATIONHOME is the inundation directory, not the data directory.
-home += sep +'data'
-
 #Derive subdirectories and filenames
-time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
-outputtimedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep+time+sep
-
+local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 
 meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
 datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
 gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
-boundarydir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
+boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
 outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
-tidedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'tide_data'+sep
+outputtimedir = outputdir + local_time + sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
 
-gauge_filename = gaugedir + 'gauge_location_broome.csv'
-buildings_filename = gaugedir + 'broome_res.csv'
-buildings_filename_out = 'broome_res_modified.csv'
-buildings_filename_damage_out = 'broome_res_modified_damage.csv'
-community_filename = gaugedir + 'CHINS_v2.csv'
-community_scenario = gaugedir + 'community_pt_hedland.csv'
-tidal_filename = tidedir + 'pt_hedland_tide.txt'
+gauge_filename = gaugedir + 'broome_gauges.csv'
 
-meshname = meshdir + basename
-onshore_dem_name = datadir + onshore_name_dli
-offshore_dem_name1 = datadir + offshore_name1
-offshore_dem_name2 = datadir + offshore_name2
+codedir = getcwd()+sep                            
+codedirname = codedir + 'project.py'
+meshname = outputtimedir + 'mesh_' + basename
+
+# Necessary if using point datasets, rather than grid
+onshore_dem_name = datadir + onshore_name
+offshore_dem_name_local1 = datadir + offshore_name_dpi1
+offshore_dem_name_aho1 = datadir + offshore_name1
+
 coast_dem_name = datadir + coast_name
-combined_dem_name = datadir + 'broome_combined_elevation'
-outputname = outputtimedir + basename  #Used by post processing
 
-# for ferret2sww
-south = degminsec2decimal_degrees(-20,30,0)
-north = degminsec2decimal_degrees(-17,10,0)
-west = degminsec2decimal_degrees(117,00,0)
-east = degminsec2decimal_degrees(120,00,0)
+combined_dem_name   = datadir + 'broome_combined_elevation'
 
-# region to export (used from export_results.py)
-e_min_area = 648000
-e_max_area = 675000
-n_min_area = 7745000
-n_max_area = 7761000
+###############################
+# Domain definitions
+###############################
 
-export_region = [[e_min_area,n_min_area],[e_min_area,n_max_area],[e_max_area,n_max_area],[e_max_area,n_min_area]]
-                 
-refzone = 50 
+# bounding box for clipping MOST/URS output (much bigger than study area)
+south = degminsec2decimal_degrees(-19,0,0)
+north = degminsec2decimal_degrees(-17,15,0)
+west  = degminsec2decimal_degrees(120,0,0)
+east  = degminsec2decimal_degrees(122,30,0) 
 
-from anuga.coordinate_transforms.redfearn import redfearn
-# boundary up to 50 m contour
-lat1_50 = degminsec2decimal_degrees(-19,20,0)
-lat2_50 = degminsec2decimal_degrees(-19,30,0)
-lat3_50 = degminsec2decimal_degrees(-19,45,0)
-lon1_50 = degminsec2decimal_degrees(119,05,0)
-lon2_50 = degminsec2decimal_degrees(118,20,0)
-lon3_50 = degminsec2decimal_degrees(117,45,0)
-z, easting, northing = redfearn(lat1_50, lon1_50)
-d0_50 = [easting, northing]
-z, easting, northing = redfearn(lat2_50, lon2_50)
-d1_50 = [easting, northing]
-z, easting, northing= redfearn(lat3_50, lon3_50)
-d2_50 = [easting, northing]
+d0 = [south, west]
+d1 = [south, east]
+d2 = [north, east]
+d3 = [north, west]
+poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
+refzone = zone
 
-d4_50 = [285000, 7585000]
-d6_50 = [330000, 7605000]
-#bounding_poly50 = [p0_50, p1_50, p2_50, d6_50, d5, d4_50]
+# bounding polygon for study area
+polyAll = read_polygon(polygondir+'extent.csv')
 
-d0 = [763852.0, 7934358.0]
-d1 = [710987.0, 7925797.0]
-d2 = [658264.0, 7926314.0]
-d3 = [552686.0, 7871580.0]
-#d4 = [604415.81, 7733013.56]
-d4 = [638000.0, 7733013.56]
-#d5 = [656561.15, 7732615.11]
-d5 = [662000.0, 7732615.11]
-#d6 = [708940.32, 7750510.33]
-d6 = [690000.0, 7740510.33]
-#polyAll = [d0, d1, d2, d3, d4, d5, d6]
-#polyAll = [d0_50, d1_50, d2_50, d4, d5, d6]
-# from Hamish
-h0=[629262.17, 7747205.47]
-h1=[552686.00, 7871579.99] #d3
-h2=[658264.00, 7926314.00] #d2
-h3=[710986.99, 7925796.99] #d1
-h4=[763851.99, 7934357.99] #d0
-h5=[701485.21, 7770656.86]
-h6=[698273.75, 7762227.38]
-h7=[698194.23, 7762018.65]
-h8=[691627.41, 7744781.98]
-h9=[679220.75, 7743604.59]
-h10=[653512.59, 7740528.56]
-h11=[634777.71, 7738247.17]
-h12=[629443.86, 7746910.37]
-h13=[629396.84, 7746986.75]
-h14=[629352.32, 7747059.06]
-h15=[629276.24, 7747182.63]
-h16=[629262.17, 7747205.47] #repeat of h0
-# using Hamish's new bounding polygon
-#polyAll = [d0_50, d1_50, d2_50, h16,h15,h14,h13,h12,h11,h10,h9,h8,h7,h6,h5]
-polyAll = [d0_50, d1_50, d2_50, h16,h11,h8,h6, h5]
+# plot bounding polygon and make sure BC info surrounds it
+plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
+print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
 
-#Interior region - Broome centre
-i0 = [668000, 7757000]
-i1 = [659000, 7755000]
-i2 = [660000, 7749000]
-i3 = [667000, 7746000]
-i4 = [678000, 7751000]
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
 
-poly_broome = [i0, i1, i2, i3, i4]
+# clipping 20m onshore data set
+#eastingmin = 520000
+#eastingmax = 536000
+#northingmin = 5245000
+#northingmax = 5260000
 
-#Are there other significant features?
-j0 = [670000, 7760000]
-j1 = [633000, 7745000]
-j2 = [665000, 7743000]
-j3 = [690000, 7755000]
+###############################
+# Interior region definitions
+###############################
 
-poly_region = [j0, j1, j2, j3]
+# broome digitized polygons
+poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon.csv')
+poly_broome2 = read_polygon(polygondir+'Broome_Close2.csv')
+poly_broome3 = read_polygon(polygondir+'Broome_Coast.csv')
+poly_broome4 = read_polygon(polygondir+'Cable_Beach.csv')
+
+plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3,poly_broome4],'boundingpoly2',verbose=False)
+print 'Area of local polygon', polygon_area(poly_broome1)/1000000.0
+print 'Area of close polygon', polygon_area(poly_broome2)/1000000.0
+print 'Area of coastal polygon', polygon_area(poly_broome3)/1000000.0
+print 'Area of cable beach polygon', polygon_area(poly_broome4)/1000000.0
+
+for i in poly_broome3:
+    v = is_inside_polygon(i,poly_broome1, verbose=False)
+    if v == False: print v
Index: /anuga_work/production/broome_2006/run_broome.py
===================================================================
--- /anuga_work/production/broome_2006/run_broome.py	(revision 3787)
+++ /anuga_work/production/broome_2006/run_broome.py	(revision 3788)
@@ -6,27 +6,39 @@
 
 The scenario is defined by a triangular mesh created from project.polygon,
-the elevation data and a simulated submarine landslide.
+the elevation data and a tsunami wave generated by MOST.
 
 Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
 """
-#-------------------------------------------------------------------------------# Import necessary modules
+def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
+    from anuga.utilities.polygon import polygon_area
+    
+    # check if any of the regions fall inside one another
+    no_triangles = 0.0
+    area = polygon_area(bounding_poly)
+    for i,j in interior_regions:
+        this_area = polygon_area(i)
+        no_triangles += this_area/j
+        area -= this_area
+        print j, this_area/1000000., area/1000000.
+    no_triangles += area/remainder_res
+    return int(no_triangles/0.7)
+#-------------------------------------------------------------------------------
+# Import necessary modules
 #-------------------------------------------------------------------------------
 
 # Standard modules
-from os import sep
-from os.path import dirname, basename
+import os
 import time
-
-# Related major packages
-from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
-                            Dirichlet_boundary, Time_boundary, File_boundary
-from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
-from anuga.pyvolution.combine_pts import combine_rectangular_points_files 
-from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
 from shutil import copy
 from os import mkdir, access, F_OK
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, \
+                            Dirichlet_boundary, Time_boundary, File_boundary
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 
 from anuga.geospatial_data.geospatial_data import *
-import sys
-from anuga.pyvolution.util import Screen_Catcher
+from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
 
 # Application specific imports
@@ -41,15 +53,16 @@
 if access(project.outputtimedir,F_OK) == 0 :
     mkdir (project.outputtimedir)
-copy (dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py')
-copy (__file__, project.outputtimedir + basename(__file__))
-print 'project.outputtimedir',project.outputtimedir
-
-# normal screen output is stored in 
+copy (project.codedirname, project.outputtimedir + project.codename)
+copy (project.codedir + 'run_broome.py', project.outputtimedir + 'run_broome.py')
+print'output dir', project.outputtimedir
+
+#normal screen output is stored in 
 screen_output_name = project.outputtimedir + "screen_output.txt"
 screen_error_name = project.outputtimedir + "screen_error.txt"
 
-# used to catch screen output to file
+#used to catch screen output to file
 sys.stdout = Screen_Catcher(screen_output_name)
 sys.stderr = Screen_Catcher(screen_error_name)
+
 print 'USER:    ', project.user
 
@@ -58,38 +71,39 @@
 #
 # Convert ASC 2 DEM 2 PTS using source data and store result in source data
-# Do for coarse and fine data
-# Fine pts file to be clipped to area of interest
 #-------------------------------------------------------------------------------
 
 # filenames
+#onshore_dem_name = project.onshore_dem_name
+#coast_points = project.coast_dem_name
+#offshore_points = project.offshore_dem_name
 meshname = project.meshname+'.msh'
-source_dir = project.boundarydir
-
-# creates DEM from asc data 
-convert_dem_from_ascii2netcdf(project.onshore_dem_name, use_cache=True, verbose=True)
-
-#creates pts file from DEM
-dem2pts(project.onshore_dem_name,
+#source_dir = project.boundarydir
+
+copied_files = False
+'''
+# creates DEM from asc data
+convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+dem2pts(onshore_dem_name,
         easting_min=project.eastingmin,
         easting_max=project.eastingmax,
         northing_min=project.northingmin,
         northing_max= project.northingmax,
-        use_cache=True,
+        use_cache=True, 
         verbose=True)
 
-print 'create G1'
-G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')
-print 'create G2'
-G2 = Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')
-print 'create G3'
-G3 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
-print 'create G4'
-G4 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
-print 'add G1+G2+G3+G4'
+print 'create offshore'
+G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya')
+print 'create onshore'
+G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
+print 'create coast'
+G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
+print 'add'
 G = G1 + G2 + G3 + G4
-print 'export G'
+print 'export points'
 G.export_points_file(project.combined_dem_name + '.pts')
-
-#-------------------------------------------------------------------------------                                 
+'''
+#----------------------------------------------------------------------------
 # Create the triangular mesh based on overall clipping polygon with a tagged
 # boundary and interior regions defined in project.py along with
@@ -98,40 +112,37 @@
 
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
-
-region_res = 500000
+remainder_res = 100000
+local_res = 15000
+broome_res = 2500
 coast_res = 500
-broome_res = 5000
-interior_regions = [[project.poly_broome, broome_res],
-                    [project.poly_region, region_res]]
+beach_res = 500
+interior_regions = [[project.poly_broome1, local_res],
+                    [project.poly_broome2, broome_res],
+                    [project.poly_broome3, coast_res],
+                    [project.poly_broome4, beach_res]]
 
 print 'number of interior regions', len(interior_regions)
-
-from anuga.utilities.polygon import plot_polygons
-if sys.platform == 'win32':
-    #figname = project.outputtimedir + 'pt_hedland_polys'
-    figname = 'broome_polys_test'
-    plot_polygons([project.polyAll,project.poly_broome,project.poly_region],
-              figname,
-              verbose = True)    
-
-print 'start create mesh from regions'
+print 'estimate of number of triangles', \
+      number_mesh_triangles(interior_regions, project.polyAll, remainder_res)
+
+import sys; sys.exit()
+
 from caching import cache
 _ = cache(create_mesh_from_regions,
           project.polyAll,
-          {'boundary_tags': {'topright': [0], 'topleft': [1],
-                             'left': [2], 'bottom0': [3],
-                             'bottom1': [4], 'bottom2': [5],
-			     'bottom3': [6], 'right': [7]},
-           'maximum_triangle_area': 250000,
-           'filename': meshname,           
+           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
+                              'e3': [3], 'e4':[4], 'e5': [5],
+			      'e6': [6]},
+           'maximum_triangle_area': remainder_res,
+           'filename': meshname,
            'interior_regions': interior_regions},
-          verbose = True, evaluate=True)
-
+          verbose = True, evaluate=False)
+print 'created mesh'
+import sys; sys.exit()
 #-------------------------------------------------------------------------------                                 
 # Setup computational domain
 #-------------------------------------------------------------------------------                                 
-domain = Domain(meshname, use_cache = False, verbose = True)
-
-print domain.statistics()
+domain = Domain(meshname, use_cache = True, verbose = True)
+
 print 'Number of triangles = ', len(domain)
 print 'The extent is ', domain.get_extent()
@@ -141,4 +152,5 @@
 domain.set_datadir(project.outputtimedir)
 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+domain.set_minimum_storable_height(0.01)
 
 #-------------------------------------------------------------------------------                                 
@@ -146,16 +158,9 @@
 #-------------------------------------------------------------------------------
 
-tide = 0.
-#high
-#tide = 
-#low
-#tide = 
-
+tide = 0.0
 domain.set_quantity('stage', tide)
 domain.set_quantity('friction', 0.0) 
-print 'hi and file',project.combined_dem_name + '.pts'
-
-domain.set_quantity('elevation', 
-                    filename = project.combined_dem_name + '.pts',
+domain.set_quantity('elevation', 0.0,
+                    #filename = project.combined_dem_name + '.pts',
                     use_cache = True,
                     verbose = True,
@@ -164,20 +169,20 @@
 
 #-------------------------------------------------------------------------------                                 
-# Setup boundary conditions 
-#-------------------------------------------------------------------------------
-print 'start ferret2sww'
-# skipped as results in file SU-AU_clipped is correct for all WA
-
-from anuga.pyvolution.data_manager import ferret2sww
+# Setup boundary conditions
+#-------------------------------------------------------------------------------
+'''
+print 'start urs2sww'
+print '', project.boundary_basename
+from anuga.shallow_water.data_manager import urs2sww
 
 south = project.south 
 north = project.north
-west = project.west
-east = project.east
+west  = project.west
+east  = project.east
 
 #note only need to do when an SWW file for the MOST boundary doesn't exist
-cache(ferret2sww,
+cache(urs2sww,
       (source_dir + project.boundary_basename,
-       source_dir + project.boundary_basename+'_'+project.basename), 
+       source_dir + project.boundary_basename), 
       {'verbose': True,
        'minlat': south,
@@ -185,23 +190,28 @@
        'minlon': west,
        'maxlon': east,
-#       'origin': project.mesh_origin,
-       'origin': domain.geo_reference.get_origin(),
+       #'origin': domain.geo_reference.get_origin(),
        'mean_stage': tide,
        'zscale': 1,                 #Enhance tsunami
        'fail_on_NaN': False,
        'inverted_bathymetry': True},
-       evaluate = True,
+      #evaluate = True,
        verbose = True,
-      dependencies = source_dir + project.boundary_basename + '.sww')
-
+       dependencies = source_dir + project.boundary_basename + '.sww')
+
+'''
 print 'Available boundary tags', domain.get_boundary_tags()
 
-Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
-                    domain, verbose = True)
+#Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
+#                    domain, verbose = True)
 Br = Reflective_boundary(domain)
 Bd = Dirichlet_boundary([tide,0,0])
-domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left':  Bd, 'bottom0': Bd,
-                      'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, 
-			'right': Bd})
+
+# 7 min square wave starting at 1 min, 6m high
+Bw = Time_boundary(domain = domain,
+                   f=lambda t: [(60<t<480)*10, 0, 0])
+
+domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
+		      'e5': Bd,  'e6': Bd} )
+
 
 #-------------------------------------------------------------------------------                                 
@@ -211,28 +221,18 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 240, finaltime = 10800): 
+for t in domain.evolve(yieldstep = 240, finaltime = 6800): 
     domain.write_time()
-    domain.write_boundary_statistics(tags = 'topright')     
-
-for t in domain.evolve(yieldstep = 120, finaltime = 16200
+    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 = 'topright')     
-
-for t in domain.evolve(yieldstep = 60, finaltime = 21600
+    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 = 'topright')     
+    domain.write_boundary_statistics(tags = 'e14') 
     
-for t in domain.evolve(yieldstep = 120, finaltime = 27000
-                       ,skip_initial_step = True): 
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'topright')     
-
-for t in domain.evolve(yieldstep = 240, finaltime = 36000
-                       ,skip_initial_step = True): 
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'topright')   
-  
 print 'That took %.2f seconds' %(time.time()-t0)
 
Index: /anuga_work/production/onslow_2006/project.py
===================================================================
--- /anuga_work/production/onslow_2006/project.py	(revision 3787)
+++ /anuga_work/production/onslow_2006/project.py	(revision 3788)
@@ -5,5 +5,5 @@
 from os import sep, environ, getenv, getcwd
 from os.path import expanduser
-#from anuga.utilities.polygon import read_polygon
+from anuga.utilities.polygon import read_polygon, polygon_area
 import sys
 
@@ -163,5 +163,5 @@
 
 polyAll = [d0, d1, d2, d3, d4, d5, d6]
-
+print 'bounding polygon area', polygon_area(polyAll)/1000000.0
 polygons = [polyAll, export_region]
 figname = 'checking.png'
@@ -208,5 +208,5 @@
 
 poly_onslow = [i0, i1, i2, i3, i4, i5]
-
+print 'onslow polygon area', polygon_area(poly_onslow)/1000000.0
 #Thevenard Island
 j0 = [294000, 7629000]
@@ -216,5 +216,5 @@
 
 poly_thevenard = [j0, j1, j2, j3]
-
+print 'thevenard polygon area', polygon_area(poly_thevenard)/1000000.0
 '''
 # Direction Is
@@ -237,5 +237,5 @@
 #poly_coast = [l0, l1, l2, l3, l4, l5]
 poly_coast = [l0, l1, l2, l3, l4]
-
+print 'coast polygon area', polygon_area(poly_coast)/1000000.0
 #general coast and local area to onslow region
 m0 = [270000, 7581000]
@@ -247,2 +247,3 @@
 
 poly_region = [m0, m1, m2, m3, m4, m5]
+print 'region polygon area', polygon_area(poly_region)/1000000.0
Index: /anuga_work/production/onslow_2006/run_onslow.py
===================================================================
--- /anuga_work/production/onslow_2006/run_onslow.py	(revision 3787)
+++ /anuga_work/production/onslow_2006/run_onslow.py	(revision 3788)
@@ -23,11 +23,10 @@
 
 # Related major packages
-from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
+from anuga.shallow_water import Domain, Reflective_boundary, \
                             Dirichlet_boundary, Time_boundary, File_boundary
-from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
-from anuga.pyvolution.combine_pts import combine_rectangular_points_files 
-from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 
 from anuga.geospatial_data.geospatial_data import *
-from anuga.pyvolution.util import Screen_Catcher
+from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
 
 # Application specific imports
@@ -144,7 +143,9 @@
 #-------------------------------------------------------------------------------                                 
 
-domain = pmesh_to_domain_instance(meshname, Domain,
-                                  use_cache = False,
-                                  verbose = True)
+#domain = pmesh_to_domain_instance(meshname, Domain,
+#                                  use_cache = False,
+#                                  verbose = True)
+
+domain = Domain(meshname, use_cache = False, verbose = True)
 
 print 'Number of triangles = ', len(domain)
Index: /anuga_work/production/pt_hedland_2006/project.py
===================================================================
--- /anuga_work/production/pt_hedland_2006/project.py	(revision 3787)
+++ /anuga_work/production/pt_hedland_2006/project.py	(revision 3788)
@@ -7,5 +7,5 @@
 #from anuga.utilities.polygon import read_polygon
 import sys
-from anuga.pmesh.create_mesh import convert_from_latlon_to_utm
+from anuga.coordinate_transforms.redfearn import convert_points_from_latlon_to_utm
 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
 from time import localtime, strftime
Index: /anuga_work/production/pt_hedland_2006/run_pt_hedland.py
===================================================================
--- /anuga_work/production/pt_hedland_2006/run_pt_hedland.py	(revision 3787)
+++ /anuga_work/production/pt_hedland_2006/run_pt_hedland.py	(revision 3788)
@@ -19,15 +19,14 @@
 
 # Related major packages
-from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
+from anuga.shallow_water import Domain, Reflective_boundary, \
                             Dirichlet_boundary, Time_boundary, File_boundary
-from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, \
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, \
      dem2pts
-from anuga.pyvolution.combine_pts import combine_rectangular_points_files 
-from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
+from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 
 from shutil import copy
 from os import mkdir, access, F_OK
 from anuga.geospatial_data.geospatial_data import *
 import sys
-from anuga.pyvolution.util import Screen_Catcher
+from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
 
 # Application specific imports
Index: /anuga_work/production/sydney_2006/project.py
===================================================================
--- /anuga_work/production/sydney_2006/project.py	(revision 3787)
+++ /anuga_work/production/sydney_2006/project.py	(revision 3788)
@@ -8,5 +8,5 @@
 import sys
 
-from anuga.pmesh.create_mesh import convert_from_latlon_to_utm
+from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_from_latlon_to_utm
                 
 
