Index: /anuga_work/production/broome_2006/locate_indigeneous.py
===================================================================
--- /anuga_work/production/broome_2006/locate_indigeneous.py	(revision 3907)
+++ /anuga_work/production/broome_2006/locate_indigeneous.py	(revision 3908)
@@ -39,7 +39,3 @@
 zone = 51
 
-print 'extent', extent
-
 community_in_extent(filename,extent,zone,fid_out)
-
-
Index: /anuga_work/production/broome_2006/project.py
===================================================================
--- /anuga_work/production/broome_2006/project.py	(revision 3907)
+++ /anuga_work/production/broome_2006/project.py	(revision 3908)
@@ -74,5 +74,6 @@
 
 gauge_filename = gaugedir + 'broome_gauges.csv'
-
+community_filename = gaugedir + 'CHINS_v2.csv'
+community_broome = gaugedir + 'community_broome.csv'
 codedir = getcwd()+sep                            
 codedirname = codedir + 'project.py'
Index: /anuga_work/production/busselton_2006/locate_indigeneous.py
===================================================================
--- /anuga_work/production/busselton_2006/locate_indigeneous.py	(revision 3908)
+++ /anuga_work/production/busselton_2006/locate_indigeneous.py	(revision 3908)
@@ -0,0 +1,41 @@
+#
+# script to determine which communities are located in study region
+# 31 October 2006
+#
+
+def community_in_extent(filename,extent,zone,fid_out):
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+
+    line1 = lines[0]
+    line11 = line1.split(',')
+    for i in range(len(line11)):
+        if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i
+        if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i
+        if line11[i].strip('\n').strip(' ') == 'COMMUNITY': name_index = i
+
+    from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
+    from anuga.utilities.polygon import plot_polygons, is_inside_polygon
+    
+    for line in lines[1:]:
+        fields = line.split(',')
+        point, zone1 = convert_from_latlon_to_utm(latitudes=[fields[lat_index]],\
+                                                 longitudes=[fields[lon_index]])
+        if zone1 == zone:
+            v = is_inside_polygon(point,extent,verbose=False)
+            if v == True:
+                p = point[0]
+                s = '%s, %.2f, %.2f\n' %(fields[name_index], p[0], p[1])
+                fid_out.write(s)
+    return
+
+import project
+
+filename = project.community_filename
+filename_out = project.community_busselton
+fid_out = open(filename_out, 'w')
+extent = project.polyAll
+zone = 50
+
+community_in_extent(filename,extent,zone,fid_out)
Index: /anuga_work/production/busselton_2006/project.py
===================================================================
--- /anuga_work/production/busselton_2006/project.py	(revision 3908)
+++ /anuga_work/production/busselton_2006/project.py	(revision 3908)
@@ -0,0 +1,174 @@
+# -*- coding: cp1252 -*-
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
+import sys
+from time import localtime, strftime, gmtime
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+
+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
+state = 'western_australia'
+scenario_dir_name = 'busselton_tsunami_scenario_2006'
+
+# onshore data provided by WA DLI
+onshore_name =  ''# original
+
+# AHO + DPI data
+offshore_name1 = 'XY100011610'
+offshore_name2 = 'XY100011611'
+offshore_name3 = 'XY100011613'
+offshore_name4 = 'XY100011614'
+offshore_name5 = 'XY100011616'
+offshore_name6 = 'XY100011617'
+offshore_name7 = 'XY100011618'
+offshore_name8 = 'XY100011621'
+offshore_name9 = 'XY100011623'
+offshore_name10 = 'XY100011745'
+offshore_name11 = 'XY100011746'
+offshore_name12 = 'XY100017530'
+offshore_name13 = 'XY100017532'
+offshore_name14 = 'XY100017538'
+offshore_name15 = 'XY100017540'
+offshore_name16 = 'XYBR66'
+offshore_name17 = 'XYBR70'
+offshore_name18 = 'XYBR80'
+offshore_name19 = 'XYBR88'
+offshore_name20 = 'XYBR93'
+offshore_name21 = 'XYBR0110'
+offshore_name22 = 'XYWADPI'
+
+# developed by NM&I
+coast_name = 'coastline'
+
+boundary_basename = 'SU-AU' # Mw ?
+
+#swollen/ all data output
+basename = 'source'
+codename = 'project.py'
+
+#Derive subdirectories and filenames
+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+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
+outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
+outputtimedir = outputdir + local_time + sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+
+gauge_filename = gaugedir + 'busselton_gauges.csv'
+community_filename = gaugedir + 'CHINS_v2.csv'
+community_broome = gaugedir + 'community_busselton.csv'
+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_name1 = datadir + offshore_name1
+offshore_dem_name2 = datadir + offshore_name2
+offshore_dem_name3 = datadir + offshore_name3
+offshore_dem_name4 = datadir + offshore_name4
+offshore_dem_name5 = datadir + offshore_name5
+offshore_dem_name6 = datadir + offshore_name6
+offshore_dem_name7 = datadir + offshore_name7
+offshore_dem_name8 = datadir + offshore_name8
+offshore_dem_name9 = datadir + offshore_name9
+offshore_dem_name10 = datadir + offshore_name10
+offshore_dem_name11 = datadir + offshore_name11
+offshore_dem_name12 = datadir + offshore_name12
+offshore_dem_name13 = datadir + offshore_name13
+offshore_dem_name14 = datadir + offshore_name14
+offshore_dem_name15 = datadir + offshore_name15
+offshore_dem_name16 = datadir + offshore_name16
+offshore_dem_name17 = datadir + offshore_name17
+offshore_dem_name18 = datadir + offshore_name18
+offshore_dem_name19 = datadir + offshore_name19
+offshore_dem_name20 = datadir + offshore_name20
+offshore_dem_name21 = datadir + offshore_name21
+offshore_dem_name22 = datadir + offshore_name22
+
+coast_dem_name = datadir + coast_name
+
+combined_dem_name   = datadir + 'broome_combined_elevation'
+
+###############################
+# Domain definitions
+###############################
+
+# 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) 
+##
+##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
+
+# bounding polygon for study area
+polyAll = read_polygon(polygondir+'extent.csv')
+
+# 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
+
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# exporting asc grid
+eastingmin = 406215.87
+eastingmax = 440208.78
+northingmin = 7983427.73
+northingmax = 8032834.52
+
+###############################
+# Interior region definitions
+###############################
+
+# broome digitized polygons
+poly_busselton1 = read_polygon(polygondir+'buss_Local_Polygon_update.csv')
+poly_busselton2 = read_polygon(polygondir+'buss_Close2_update.csv')
+poly_busselton3 = read_polygon(polygondir+'buss_Coast_update.csv')
+
+plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],'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
+
+for i in poly_broome3:
+    v = is_inside_polygon(i,poly_broome1, verbose=False)
+    if v == False: print v
+
+def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
+    from anuga.utilities.polygon import polygon_area
+    
+    # TO DO 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)
Index: /anuga_work/production/exmouth_2006/locate_indigeneous.py
===================================================================
--- /anuga_work/production/exmouth_2006/locate_indigeneous.py	(revision 3908)
+++ /anuga_work/production/exmouth_2006/locate_indigeneous.py	(revision 3908)
@@ -0,0 +1,41 @@
+#
+# script to determine which communities are located in study region
+# 31 October 2006
+#
+
+def community_in_extent(filename,extent,zone,fid_out):
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+
+    line1 = lines[0]
+    line11 = line1.split(',')
+    for i in range(len(line11)):
+        if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i
+        if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i
+        if line11[i].strip('\n').strip(' ') == 'COMMUNITY': name_index = i
+
+    from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
+    from anuga.utilities.polygon import plot_polygons, is_inside_polygon
+    
+    for line in lines[1:]:
+        fields = line.split(',')
+        point, zone1 = convert_from_latlon_to_utm(latitudes=[fields[lat_index]],\
+                                                 longitudes=[fields[lon_index]])
+        if zone1 == zone:
+            v = is_inside_polygon(point,extent,verbose=False)
+            if v == True:
+                p = point[0]
+                s = '%s, %.2f, %.2f\n' %(fields[name_index], p[0], p[1])
+                fid_out.write(s)
+    return
+
+import project
+
+filename = project.community_filename
+filename_out = project.community_broome
+fid_out = open(filename_out, 'w')
+extent = project.polyAll
+zone = 50 # or is it 49?
+
+community_in_extent(filename,extent,zone,fid_out)
Index: /anuga_work/production/exmouth_2006/project.py
===================================================================
--- /anuga_work/production/exmouth_2006/project.py	(revision 3908)
+++ /anuga_work/production/exmouth_2006/project.py	(revision 3908)
@@ -0,0 +1,177 @@
+# -*- coding: cp1252 -*-
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
+import sys
+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
+state = 'western_australia'
+scenario_dir_name = 'exmouth_tsunami_scenario_2006'
+
+# onshore data provided by WA DLI
+#onshore_name =  # original
+
+# AHO + DPI data
+offshore_name1 = 'XY100011610'
+offshore_name2 = 'XY100011611'
+offshore_name3 = 'XY100011613'
+offshore_name4 = 'XY100011614'
+offshore_name5 = 'XY100011616'
+offshore_name6 = 'XY100011617'
+offshore_name7 = 'XY100011618'
+offshore_name8 = 'XY100011621'
+offshore_name9 = 'XY100011623'
+offshore_name10 = 'XY100011745'
+offshore_name11 = 'XY100011746'
+offshore_name12 = 'XY100017530'
+offshore_name13 = 'XY100017532'
+offshore_name14 = 'XY100017538'
+offshore_name15 = 'XY100017540'
+offshore_name16 = 'XYBR66'
+offshore_name17 = 'XYBR70'
+offshore_name18 = 'XYBR80'
+offshore_name19 = 'XYBR88'
+offshore_name20 = 'XYBR93'
+offshore_name21 = 'XYBR0110'
+offshore_name22 = 'XYWADPI'
+
+# developed by NM&I
+coast_name = 'coastline'
+
+boundary_basename = 'SU-AU' # Mw ?
+
+#swollen/ all data output
+basename = 'source'
+codename = 'project.py'
+
+#Derive subdirectories and filenames
+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+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
+outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
+outputtimedir = outputdir + local_time + sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+
+gauge_filename = gaugedir + 'exmouth_gauges.csv'
+community_filename = gaugedir + 'CHINS_v2.csv'
+community_exmouth = gaugedir + 'community_exmouth.csv'
+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_name1 = datadir + offshore_name1
+offshore_dem_name2 = datadir + offshore_name2
+offshore_dem_name3 = datadir + offshore_name3
+offshore_dem_name4 = datadir + offshore_name4
+offshore_dem_name5 = datadir + offshore_name5
+offshore_dem_name6 = datadir + offshore_name6
+offshore_dem_name7 = datadir + offshore_name7
+offshore_dem_name8 = datadir + offshore_name8
+offshore_dem_name9 = datadir + offshore_name9
+offshore_dem_name10 = datadir + offshore_name10
+offshore_dem_name11 = datadir + offshore_name11
+offshore_dem_name12 = datadir + offshore_name12
+offshore_dem_name13 = datadir + offshore_name13
+offshore_dem_name14 = datadir + offshore_name14
+offshore_dem_name15 = datadir + offshore_name15
+offshore_dem_name16 = datadir + offshore_name16
+offshore_dem_name17 = datadir + offshore_name17
+offshore_dem_name18 = datadir + offshore_name18
+offshore_dem_name19 = datadir + offshore_name19
+offshore_dem_name20 = datadir + offshore_name20
+offshore_dem_name21 = datadir + offshore_name21
+offshore_dem_name22 = datadir + offshore_name22
+
+coast_dem_name = datadir + coast_name
+
+combined_dem_name   = datadir + 'exmouth_combined_elevation'
+
+###############################
+# Domain definitions
+###############################
+
+# 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) 
+##
+##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
+
+# bounding polygon for study area
+polyAll = read_polygon(polygondir+'extent.csv')
+
+# 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
+
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# exporting asc grid
+eastingmin = 406215.87
+eastingmax = 440208.78
+northingmin = 7983427.73
+northingmax = 8032834.52
+
+###############################
+# Interior region definitions
+###############################
+
+# exmouth digitized polygons
+poly_exmouth1 = read_polygon(polygondir+'exmouth_Local_Polygon_update.csv')
+poly_exmouth2 = read_polygon(polygondir+'exmouth_Close2_update.csv')
+poly_exmouth3 = read_polygon(polygondir+'exmouth_Coast_update.csv')
+#poly_exmouth4 = read_polygon(polygondir+'Cable_Beach_revised.csv')
+
+plot_polygons([polyAll,poly_exmouth1,poly_exmouth2,poly_exmouth3],'boundingpoly2',verbose=False)
+print 'Area of local polygon', polygon_area(poly_exmouth1)/1000000.0
+print 'Area of close polygon', polygon_area(poly_exmouth2)/1000000.0
+print 'Area of coastal polygon', polygon_area(poly_exmouth3)/1000000.0
+#print 'Area of cable beach polygon', polygon_area(poly_exmouth4)/1000000.0
+
+for i in poly_exmouth3:
+    v = is_inside_polygon(i,poly_exmouth1, verbose=False)
+    if v == False: print v
+
+def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
+    from anuga.utilities.polygon import polygon_area
+    
+    # TO DO 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)
Index: /anuga_work/production/perth_2006/locate_indigeneous.py
===================================================================
--- /anuga_work/production/perth_2006/locate_indigeneous.py	(revision 3908)
+++ /anuga_work/production/perth_2006/locate_indigeneous.py	(revision 3908)
@@ -0,0 +1,41 @@
+#
+# script to determine which communities are located in study region
+# 31 October 2006
+#
+
+def community_in_extent(filename,extent,zone,fid_out):
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+
+    line1 = lines[0]
+    line11 = line1.split(',')
+    for i in range(len(line11)):
+        if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i
+        if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i
+        if line11[i].strip('\n').strip(' ') == 'COMMUNITY': name_index = i
+
+    from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
+    from anuga.utilities.polygon import plot_polygons, is_inside_polygon
+    
+    for line in lines[1:]:
+        fields = line.split(',')
+        point, zone1 = convert_from_latlon_to_utm(latitudes=[fields[lat_index]],\
+                                                 longitudes=[fields[lon_index]])
+        if zone1 == zone:
+            v = is_inside_polygon(point,extent,verbose=False)
+            if v == True:
+                p = point[0]
+                s = '%s, %.2f, %.2f\n' %(fields[name_index], p[0], p[1])
+                fid_out.write(s)
+    return
+
+import project
+
+filename = project.community_filename
+filename_out = project.community_perth
+fid_out = open(filename_out, 'w')
+extent = project.polyAll
+zone = 50
+
+community_in_extent(filename,extent,zone,fid_out)
Index: /anuga_work/production/perth_2006/project.py
===================================================================
--- /anuga_work/production/perth_2006/project.py	(revision 3908)
+++ /anuga_work/production/perth_2006/project.py	(revision 3908)
@@ -0,0 +1,177 @@
+# -*- coding: cp1252 -*-
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
+import sys
+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
+state = 'western_australia'
+scenario_dir_name = 'perth_tsunami_scenario_2006'
+
+# onshore data provided by WA DLI
+onshore_name = '' # original
+
+# AHO + DPI data
+offshore_name1 = 'XY100011610'
+offshore_name2 = 'XY100011611'
+offshore_name3 = 'XY100011613'
+offshore_name4 = 'XY100011614'
+offshore_name5 = 'XY100011616'
+offshore_name6 = 'XY100011617'
+offshore_name7 = 'XY100011618'
+offshore_name8 = 'XY100011621'
+offshore_name9 = 'XY100011623'
+offshore_name10 = 'XY100011745'
+offshore_name11 = 'XY100011746'
+offshore_name12 = 'XY100017530'
+offshore_name13 = 'XY100017532'
+offshore_name14 = 'XY100017538'
+offshore_name15 = 'XY100017540'
+offshore_name16 = 'XYBR66'
+offshore_name17 = 'XYBR70'
+offshore_name18 = 'XYBR80'
+offshore_name19 = 'XYBR88'
+offshore_name20 = 'XYBR93'
+offshore_name21 = 'XYBR0110'
+offshore_name22 = 'XYWADPI'
+
+# developed by NM&I
+coast_name = 'coastline'
+
+boundary_basename = 'SU-AU' # Mw ?
+
+#swollen/ all data output
+basename = 'source'
+codename = 'project.py'
+
+#Derive subdirectories and filenames
+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+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
+outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
+outputtimedir = outputdir + local_time + sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+
+gauge_filename = gaugedir + 'perth_gauges.csv'
+community_filename = gaugedir + 'CHINS_v2.csv'
+community_perth = gaugedir + 'community_perth.csv'
+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_name1 = datadir + offshore_name1
+offshore_dem_name2 = datadir + offshore_name2
+offshore_dem_name3 = datadir + offshore_name3
+offshore_dem_name4 = datadir + offshore_name4
+offshore_dem_name5 = datadir + offshore_name5
+offshore_dem_name6 = datadir + offshore_name6
+offshore_dem_name7 = datadir + offshore_name7
+offshore_dem_name8 = datadir + offshore_name8
+offshore_dem_name9 = datadir + offshore_name9
+offshore_dem_name10 = datadir + offshore_name10
+offshore_dem_name11 = datadir + offshore_name11
+offshore_dem_name12 = datadir + offshore_name12
+offshore_dem_name13 = datadir + offshore_name13
+offshore_dem_name14 = datadir + offshore_name14
+offshore_dem_name15 = datadir + offshore_name15
+offshore_dem_name16 = datadir + offshore_name16
+offshore_dem_name17 = datadir + offshore_name17
+offshore_dem_name18 = datadir + offshore_name18
+offshore_dem_name19 = datadir + offshore_name19
+offshore_dem_name20 = datadir + offshore_name20
+offshore_dem_name21 = datadir + offshore_name21
+offshore_dem_name22 = datadir + offshore_name22
+
+coast_dem_name = datadir + coast_name
+
+combined_dem_name   = datadir + 'perth_combined_elevation'
+
+###############################
+# Domain definitions
+###############################
+
+# 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) 
+##
+##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
+
+# bounding polygon for study area
+polyAll = read_polygon(polygondir+'extent.csv')
+
+# 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
+
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# exporting asc grid
+eastingmin = 406215.87
+eastingmax = 440208.78
+northingmin = 7983427.73
+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
+
+for i in poly_perth3:
+    v = is_inside_polygon(i,poly_perth1, verbose=False)
+    if v == False: print v
+
+def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
+    from anuga.utilities.polygon import polygon_area
+    
+    # TO DO 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)
