Index: /anuga_work/production/dampier_2006/build_dampier.py
===================================================================
--- /anuga_work/production/dampier_2006/build_dampier.py	(revision 3827)
+++ /anuga_work/production/dampier_2006/build_dampier.py	(revision 3827)
@@ -0,0 +1,178 @@
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.output_time_dir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated submarine landslide.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+from os import sep
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
+import time
+import sys
+
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import File_boundary
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+from anuga.geospatial_data.geospatial_data import *
+
+# Application specific imports
+import project   # Definition of file names and polygons
+
+
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+
+
+# creates copy of code in output dir if dir doesn't exist
+if access(project.output_time_dir,F_OK) == 0:
+    mkdir (project.output_time_dir)
+copy (dirname(project.__file__) +sep+ project.__name__+'.py',
+      project.output_time_dir + project.__name__+'.py')     #copies project.py
+copy (__file__, project.output_time_dir + basename(__file__))
+print 'project.output_time_dir',project.output_time_dir       #copies this file
+
+
+
+
+boundary_dir_name = project.boundary_dir_name
+
+
+from anuga.shallow_water.data_manager import urs2sww
+
+urs2sww(boundary_dir_name,
+#        minlat=project.north_boundary, maxlat=project.south_boundary,
+#        minlon= project.west_boundary, maxlon=project.east_boundary,
+        verbose='true')
+        
+import sys; sys.exit()
+#-------------------------------------------------------------------------------
+# Preparation of topographic data
+#
+# 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
+#-------------------------------------------------------------------------------
+
+# topography directory filenames
+onshore_dir_name = project.onshore_dir_name
+coast_dir_name = project.coast_dir_name
+islands_dir_name = project.islands_dir_name
+offshore_dir_name = project.offshore_dir_name
+offshore_dir_name1 = project.offshore_dir_name1
+offshore_dir_name2 = project.offshore_dir_name2
+offshore_dir_name3 = project.offshore_dir_name3
+offshore_dir_name4 = project.offshore_dir_name4
+offshore_dir_name5 = project.offshore_dir_name5
+offshore_dir_name6 = project.offshore_dir_name6
+offshore_dir_name7 = project.offshore_dir_name7
+offshore_dir_name8 = project.offshore_dir_name8
+offshore_dir_name9 = project.offshore_dir_name9
+offshore_dir_name10 = project.offshore_dir_name10
+offshore_dir_name11 = project.offshore_dir_name11
+offshore_dir_name12 = project.offshore_dir_name12
+offshore_dir_name13 = project.offshore_dir_name13
+offshore_dir_name14 = project.offshore_dir_name14
+
+
+
+# files to be used
+#file_used = []
+
+# creates DEM from asc data
+convert_dem_from_ascii2netcdf(onshore_dir_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(islands_dir_name, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+dem2pts(onshore_dir_name,
+#        easting_min=project.eastingmin,
+#        easting_max=project.eastingmax,
+#        northing_min=project.northingmin,
+#        northing_max= project.northingmax,
+        use_cache=True, 
+        verbose=True)
+
+
+#creates pts file for islands DEM
+dem2pts(islands_dir_name, use_cache=True, verbose=True)
+
+
+print'create G1'
+G1 = Geospatial_data(file_name = project.onshore_dir_name + '.pts')
+G2 = Geospatial_data(file_name = project.coast_dir_name + '.xya')
+G3 = Geospatial_data(file_name = project.islands_dir_name + '.pts')
+G_off = Geospatial_data(file_name = project.offshore_dir_name + '.xya')
+G_off1 = Geospatial_data(file_name = project.offshore_dir_name1 + '.xya')
+G_off2 = Geospatial_data(file_name = project.offshore_dir_name2 + '.xya')
+G_off3 = Geospatial_data(file_name = project.offshore_dir_name3 + '.xya')
+G_off4 = Geospatial_data(file_name = project.offshore_dir_name4 + '.xya')
+G_off5 = Geospatial_data(file_name = project.offshore_dir_name5 + '.xya')
+G_off6 = Geospatial_data(file_name = project.offshore_dir_name6 + '.xya')
+G_off7 = Geospatial_data(file_name = project.offshore_dir_name7 + '.xya')
+G_off8 = Geospatial_data(file_name = project.offshore_dir_name8 + '.xya')
+G_off9 = Geospatial_data(file_name = project.offshore_dir_name9 + '.xya')
+G_off10 = Geospatial_data(file_name = project.offshore_dir_name10 + '.xya')
+G_off11 = Geospatial_data(file_name = project.offshore_dir_name11 + '.xya')
+G_off12 = Geospatial_data(file_name = project.offshore_dir_name12 + '.xya')
+G_off13 = Geospatial_data(file_name = project.offshore_dir_name13 + '.xya')
+G_off14 = Geospatial_data(file_name = project.offshore_dir_name14 + '.xya')
+
+print'add G1+G2+G3+all offshore data'
+G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 \
+    + G_off6 + G_off7 + G_off8 + G_off9 + G_off10 + G_off11 + G_off12 \
+    + G_off13 + G_off14
+
+G.clip(project.bounding_polygon)
+#FIXME: add a clip function to pts
+
+print'export G'
+G.export_points_file(project.combined_dir_name + '.pts')
+import sys; sys.exit()
+
+
+
+#-------------------------------------------------------------------------
+# Convert URS to SWW file for boundary conditions 
+#-------------------------------------------------------------------------
+
+# filenames
+meshname = project.mesh_name+'.msh'
+source_dir = project.boundarydir
+boundary_file = project.boundaryname
+
+print 'Available boundary tags', domain.get_boundary_tags()
+
+from anuga.shallow_water.data_manager import urs2sww
+
+urs2sww(boundary_file,
+        minlat=project.north_boundary, maxlat=project.south_boundary,
+        minlon= project.west_boundary, maxlon=project.east_boundary,
+        verbose='true')
+
+
+
+
+
+
Index: /anuga_work/production/dampier_2006/export_results.py
===================================================================
--- /anuga_work/production/dampier_2006/export_results.py	(revision 3827)
+++ /anuga_work/production/dampier_2006/export_results.py	(revision 3827)
@@ -0,0 +1,76 @@
+import project, os
+import sys
+
+from anuga.shallow_water.data_manager import sww2dem
+#from anuga.pyvolution.ermapper_grids import read_ermapper_grid
+from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
+from os import sep
+
+#time_dir = '20060704_063005' #HAT
+#time_dir = '20060706_235246' #LAT
+time_dir = '20060704_063234' #MSL
+#time_dir = '20060515_001733' #DTED data
+directory = project.boundarydir
+#name = directory + time_dir +sep + 'source'
+# to develop grid of MOST output
+#directory = project.boundarydir
+name = directory + project.boundary_basename
+
+#normal screen output is stored in 'outname'
+#screen_output_name = directory + time_dir + sep + "export_output.txt"
+#screen_error_name = directory + time_dir + sep + "export_error.txt"
+# for MOST output
+#screen_output_name = directory  +  "export_output.txt"
+#screen_error_name = directory +  "export_error.txt"
+
+#used to catch screen output to file
+#sys.stdout = Screen_Catcher(screen_output_name)
+#sys.stderr = Screen_Catcher(screen_error_name)
+
+print 'output dir:', name
+#print 'Which variable do you want to export?'
+#which_var = int(raw_input('Stage = 0, Absolute Momentum = 1, Depth = 2, Speed = 3  '  ))
+which_var = 4
+#sys.stderr.write(sys.stdout.data)
+if which_var == 0:  # Stage
+    outname = name + '_stage'
+    quantityname = 'stage'
+
+if which_var == 1:  # Absolute Momentum
+    outname = name + '_momentum'
+    quantityname = '(xmomentum**2 + ymomentum**2)**0.5'  #Absolute momentum
+
+if which_var == 2:  # Depth
+    outname = name + '_depth'
+    quantityname = 'stage-elevation'  #Depth
+
+if which_var == 3:  # Speed
+    outname = name + '_speed'
+    #quantityname = '((xmomentum/(stage-elevation))**2 + (ymomentum/(stage-elevation))**2)**0.5'  #Speed
+    #quantityname = 'xmomentum/(stage-elevation)'  #Speed
+    quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-30)'  #Speed
+
+if which_var == 4:  # Elevation
+    outname = name + '_elevation'
+    quantityname = 'elevation'  #Elevation
+
+print 'start sww2dem'
+#sys.stderr.write(sys.stdout.data)
+sww2dem(name, basename_out = outname,
+            quantity = quantityname,
+            cellsize = 20,       # Trevor would like this at 25
+            # define region for viz purposes
+            easting_min = project.e_min_area,
+            easting_max = project.e_max_area,
+            northing_min = project.n_min_area,
+            northing_max = project.n_max_area,        
+            reduction = max, #this is because we want max quantityname
+            verbose = True,
+            format = 'asc')
+
+#sys.stderr.write(sys.stdout.data)
+
+#Check
+
+#data = read_ermapper_grid(name)
+#print 'Values from %s are in [%f, %f]' %(name, min(data.flat), max(data.flat))
Index: /anuga_work/production/dampier_2006/project_new.py
===================================================================
--- /anuga_work/production/dampier_2006/project_new.py	(revision 3827)
+++ /anuga_work/production/dampier_2006/project_new.py	(revision 3827)
@@ -0,0 +1,277 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+import sys
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser, basename
+
+from anuga.coordinate_transforms.redfearn import\
+     degminsec2decimal_degrees,\
+     convert_from_latlon_to_utm
+
+from time import localtime, strftime
+from anuga.geospatial_data.geospatial_data import *
+
+# file and system info
+#---------------------------------
+codename = 'project.py'
+
+home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
+
+if sys.platform == 'win32':
+    user = getenv('USERPROFILE')
+else:
+    user = getenv('LOGNAME')
+
+# INUNDATIONHOME is the inundation directory, not the data directory.
+home += sep +'data'
+#----------------------------------
+# Location and naming of scenario data
+#----------------------------------
+state = 'western_australia'
+scenario_name = 'dampier_tsunami'
+scenario_datas_name = 'dampier_tsunami_scenario_2006'  #name of the directory where the data is stored
+#scenario_datas_name = 'karratha_tsunami_scenario_2005' # Tmp location
+
+mesh_name = 'elevation50m'
+boundary_name = 'dampier'
+
+# topography file names
+onshore_name = 'dli_no_islands.asc'
+coast_name = 'DTED_05_Contour.xya'
+islands_name = 'dted_islands.asc'
+offshore_name = 'XY100003902.xya'
+offshore_name1 = 'XY100003903.xya'
+offshore_name2 = 'XY100003951.xya'
+offshore_name3 = 'XY100006321.xya'
+offshore_name4 = 'XY100011756.xya'
+offshore_name5 = 'XY100014243.xya'
+offshore_name6 = 'XY100014244.xya'
+offshore_name7 = 'XY100021081.xya'
+offshore_name8 = 'XY100021082.xya'
+offshore_name9 = 'XY100021083.xya'
+offshore_name10 = 'XY100021085.xya'
+offshore_name11 = 'XY100021086.xya'
+offshore_name12 = 'XY100026309.xya'
+offshore_name13 = 'XY100026338.xya'
+offshore_name14 = 'XYDM83.xya'
+
+combined_name ='dampier_combined_elevation'
+
+
+
+#Derive subdirectories and filenames
+
+meshes_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'meshes'+sep
+topographies_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'topographies'+sep
+gauges_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep
+polygons_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'polygons'+sep
+boundaries_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+'test'+sep
+#outputdir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'output'+sep
+tide_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'tide_data'+sep
+
+time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
+
+output_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'output'+sep+time+sep
+topographies_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'topographies'+sep+time+sep
+boundaries_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+time+sep
+
+
+
+#gauge_filename = gaugedir + 'gauge_location_broome.csv'
+
+onshore_dir_name = topographies_dir + onshore_name
+coast_dir_name = topographies_dir + coast_name
+islands_dir_name = topographies_dir + islands_name
+offshore_dir_name = topographies_dir + offshore_name
+
+combined_dir_name = topographies_dir + combined_name
+
+mesh_dir_name = mesh_dir + scenario_name
+output_dir_name = output_time_dir + scenario_name  #Used by post processing
+boundary_dir_name = boundary_dir + boundary_name  #Used by post processing
+
+
+
+
+# Regions
+
+refzone = 50 
+south = degminsec2decimal_degrees(-20,55,0)
+north = degminsec2decimal_degrees(-19,15,0)
+west = degminsec2decimal_degrees(116,17,0)
+east = degminsec2decimal_degrees(117,10,0)
+
+p0 = [south, degminsec2decimal_degrees(116,32,0)]
+p1 = [south, west]
+p2 = [degminsec2decimal_degrees(-20,23,0), west]
+p3 = [north, degminsec2decimal_degrees(116,45,0)]
+p4 = [north, degminsec2decimal_degrees(117,0,0)]
+p5 = [p2[0], degminsec2decimal_degrees(117,8,0)]
+p6 = [degminsec2decimal_degrees(-20,30,0), east]
+p7 = [degminsec2decimal_degrees(-20,38,0), east]
+p8 = [south, east]
+
+
+bounding_polygon, zone =\
+                  convert_from_latlon_to_utm([p0, p1, p2, p3, p4,
+                                              p5, p6, p7, p8])
+refzone = zone
+
+
+#Interior regions
+
+# CIPMA point of interest
+cipma_latitude = -20.588456
+cipma_longitude = 116.771527
+
+
+k0 = [cipma_latitude-0.02, cipma_longitude-0.02]
+k1 = [cipma_latitude-0.02, cipma_longitude+0.02]
+k2 = [cipma_latitude+0.02, cipma_longitude+0.02]
+k3 = [cipma_latitude+0.02, cipma_longitude-0.02]
+
+cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
+assert zone == refzone
+
+e_min_area = 474000
+e_max_area = 480000
+n_min_area = 7719000
+n_max_area = 7725000
+
+#Interior regions
+karratha_south = degminsec2decimal_degrees(-20,44,0)
+karratha_north = degminsec2decimal_degrees(-20,42,0)
+karratha_west = degminsec2decimal_degrees(116,48,0)
+karratha_east = degminsec2decimal_degrees(116,53,30)
+
+k0 = [karratha_south, karratha_west]
+k1 = [karratha_south, karratha_east]
+k2 = [karratha_north, karratha_east]
+k3 = [karratha_north, karratha_west]    
+
+karratha_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
+assert zone == refzone
+
+
+#Interior regions
+dampier_south = degminsec2decimal_degrees(-20,40,0)
+dampier_north = degminsec2decimal_degrees(-20,38,10)
+dampier_west = degminsec2decimal_degrees(116,43,0)
+dampier_east = degminsec2decimal_degrees(116,45,0)
+
+d0 = [dampier_south, dampier_west]
+d1 = [dampier_south, dampier_east]
+d2 = [dampier_north, dampier_east]
+d3 = [dampier_north, dampier_west]    
+
+dampier_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
+assert zone == refzone
+
+
+#Interior regions
+refinery_south = degminsec2decimal_degrees(-20,37,50)
+refinery_north = degminsec2decimal_degrees(-20,36,0)
+refinery_west = degminsec2decimal_degrees(116,44,0)
+refinery_east = degminsec2decimal_degrees(116,46,10)
+
+d0 = [refinery_south, refinery_west]
+d1 = [refinery_south, refinery_east]
+d2 = [refinery_north, refinery_east]
+d3 = [refinery_north, refinery_west]    
+
+refinery_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
+assert zone == refzone
+
+
+#Interior region around 468899, 7715177:
+#lat (-20, 39, 44.93753), lon (116, 42, 5.09106)
+
+point_south = degminsec2decimal_degrees(-20,39,46)
+point_north = degminsec2decimal_degrees(-20,39,42)
+point_west = degminsec2decimal_degrees(116,42,0)
+point_east = degminsec2decimal_degrees(116,42,10)
+
+d0 = [point_south, point_west]
+d1 = [point_south, point_east]
+d2 = [point_north, point_east]
+d3 = [point_north, point_west]    
+
+point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
+assert zone == refzone
+
+
+#Neils areas around interesting points
+neil1_point1 = [degminsec2decimal_degrees(-20,35,34),
+                degminsec2decimal_degrees(116,45,18)]
+neil1_point2 = [degminsec2decimal_degrees(-20,36,15),
+                degminsec2decimal_degrees(116,46,18)]
+neil1_point3 = [degminsec2decimal_degrees(-20,35,9),
+                degminsec2decimal_degrees(116,47,17)]
+neil1_point4 = [degminsec2decimal_degrees(-20,34,26),
+                degminsec2decimal_degrees(116,46,17)]
+
+neil1_polygon, zone = convert_from_latlon_to_utm([neil1_point1,
+                                                  neil1_point2,
+                                                  neil1_point3,
+                                                  neil1_point4])
+assert zone == refzone
+
+
+
+
+neil2_point1 = [degminsec2decimal_degrees(-20,39,36),
+                degminsec2decimal_degrees(116,41,33)]
+neil2_point2 = [degminsec2decimal_degrees(-20,40,10),
+                degminsec2decimal_degrees(116,42,13)]
+neil2_point3 = [degminsec2decimal_degrees(-20,38,39),
+                degminsec2decimal_degrees(116,43,49)]
+neil2_point4 = [degminsec2decimal_degrees(-20,38,5),
+                degminsec2decimal_degrees(116,43,9)]
+
+neil2_polygon, zone = convert_from_latlon_to_utm([neil2_point1,
+                                                  neil2_point2,
+                                                  neil2_point3,
+                                                  neil2_point4])
+assert zone == refzone
+
+
+
+
+
+#Withnell bay
+wb_point1 = [degminsec2decimal_degrees(-20,35,34),
+                degminsec2decimal_degrees(116,45,18)]
+wb_point2 = [degminsec2decimal_degrees(-20,36,15),
+                degminsec2decimal_degrees(116,46,18)]
+wb_point3 = [degminsec2decimal_degrees(-20,35,9),
+                degminsec2decimal_degrees(116,47,17)]
+wb_point4 = [degminsec2decimal_degrees(-20,34,26),
+                degminsec2decimal_degrees(116,46,17)]
+
+wb_polygon, zone = convert_from_latlon_to_utm([wb_point1, wb_point2,
+                                               wb_point3, wb_point4])
+assert zone == refzone
+
+
+
+
+
+#Larger Withnell bay
+lwb_point1 = [degminsec2decimal_degrees(-20,35,59),
+                degminsec2decimal_degrees(116,42,00)]
+lwb_point2 = [degminsec2decimal_degrees(-20,36,50),
+                degminsec2decimal_degrees(116,46,50)]
+lwb_point3 = [degminsec2decimal_degrees(-20,34,00),
+                degminsec2decimal_degrees(116,47,39)]
+lwb_point4 = [degminsec2decimal_degrees(-20,33,00),
+                degminsec2decimal_degrees(116,42,50)]
+
+lwb_polygon, zone = convert_from_latlon_to_utm([lwb_point1, lwb_point2,
+                                                lwb_point3, lwb_point4])
+                                                      
+assert zone == refzone
+
+
+
+
Index: /anuga_work/production/dampier_2006/project_original.py
===================================================================
--- /anuga_work/production/dampier_2006/project_original.py	(revision 3827)
+++ /anuga_work/production/dampier_2006/project_original.py	(revision 3827)
@@ -0,0 +1,239 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+import sys
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser, basename
+
+from anuga.coordinate_transforms.redfearn import\
+     degminsec2decimal_degrees,\
+     convert_from_latlon_to_utm
+
+from time import localtime, strftime
+from anuga.geospatial_data.geospatial_data import *
+
+# Location and naming of scenario data
+state = 'western_australia'
+scenario_dir_name = 'dampier_tsunami_scenario_2006'
+#scenario_dir_name = 'karratha_tsunami_scenario_2005' # Tmp location
+
+basename = 'elevation50m'
+boundary_basename = 'dampier'
+
+
+# Data output
+codename = 'project.py'
+
+home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
+
+if sys.platform == 'win32':
+    user = getenv('USERPROFILE')
+else:
+    user = getenv('LOGNAME')
+
+# INUNDATIONHOME is the inundation directory, not the data directory.
+home += sep +'data'
+
+#Derive subdirectories and filenames
+print 'home',home
+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
+
+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+'urs'+sep+'test'+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
+
+#gauge_filename = gaugedir + 'gauge_location_broome.csv'
+
+meshname = meshdir + basename
+outputname = outputtimedir + basename  #Used by post processing
+boundaryname = boundarydir + boundary_basename  #Used by post processing
+
+
+
+
+# Regions
+
+refzone = 50 
+south = degminsec2decimal_degrees(-20,55,0)
+north = degminsec2decimal_degrees(-20,15,0)
+west = degminsec2decimal_degrees(116,17,0)
+east = degminsec2decimal_degrees(117,10,0)
+
+p0 = [south, degminsec2decimal_degrees(116,32,0)]
+p1 = [south, west]
+p2 = [degminsec2decimal_degrees(-20,23,0), west]
+p3 = [north, degminsec2decimal_degrees(116,45,0)]
+p4 = [north, degminsec2decimal_degrees(117,0,0)]
+p5 = [p2[0], degminsec2decimal_degrees(117,8,0)]
+p6 = [degminsec2decimal_degrees(-20,30,0), east]
+p7 = [degminsec2decimal_degrees(-20,38,0), east]
+p8 = [south, east]
+
+
+bounding_polygon, zone =\
+                  convert_from_latlon_to_utm([p0, p1, p2, p3, p4,
+                                              p5, p6, p7, p8])
+refzone = zone
+
+
+#Interior regions
+
+# CIPMA point of interest
+cipma_latitude = -20.588456
+cipma_longitude = 116.771527
+
+
+k0 = [cipma_latitude-0.02, cipma_longitude-0.02]
+k1 = [cipma_latitude-0.02, cipma_longitude+0.02]
+k2 = [cipma_latitude+0.02, cipma_longitude+0.02]
+k3 = [cipma_latitude+0.02, cipma_longitude-0.02]
+
+cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
+assert zone == refzone
+
+e_min_area = 474000
+e_max_area = 480000
+n_min_area = 7719000
+n_max_area = 7725000
+
+#Interior regions
+karratha_south = degminsec2decimal_degrees(-20,44,0)
+karratha_north = degminsec2decimal_degrees(-20,42,0)
+karratha_west = degminsec2decimal_degrees(116,48,0)
+karratha_east = degminsec2decimal_degrees(116,53,30)
+
+k0 = [karratha_south, karratha_west]
+k1 = [karratha_south, karratha_east]
+k2 = [karratha_north, karratha_east]
+k3 = [karratha_north, karratha_west]    
+
+karratha_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
+assert zone == refzone
+
+
+#Interior regions
+dampier_south = degminsec2decimal_degrees(-20,40,0)
+dampier_north = degminsec2decimal_degrees(-20,38,10)
+dampier_west = degminsec2decimal_degrees(116,43,0)
+dampier_east = degminsec2decimal_degrees(116,45,0)
+
+d0 = [dampier_south, dampier_west]
+d1 = [dampier_south, dampier_east]
+d2 = [dampier_north, dampier_east]
+d3 = [dampier_north, dampier_west]    
+
+dampier_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
+assert zone == refzone
+
+
+#Interior regions
+refinery_south = degminsec2decimal_degrees(-20,37,50)
+refinery_north = degminsec2decimal_degrees(-20,36,0)
+refinery_west = degminsec2decimal_degrees(116,44,0)
+refinery_east = degminsec2decimal_degrees(116,46,10)
+
+d0 = [refinery_south, refinery_west]
+d1 = [refinery_south, refinery_east]
+d2 = [refinery_north, refinery_east]
+d3 = [refinery_north, refinery_west]    
+
+refinery_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
+assert zone == refzone
+
+
+#Interior region around 468899, 7715177:
+#lat (-20, 39, 44.93753), lon (116, 42, 5.09106)
+
+point_south = degminsec2decimal_degrees(-20,39,46)
+point_north = degminsec2decimal_degrees(-20,39,42)
+point_west = degminsec2decimal_degrees(116,42,0)
+point_east = degminsec2decimal_degrees(116,42,10)
+
+d0 = [point_south, point_west]
+d1 = [point_south, point_east]
+d2 = [point_north, point_east]
+d3 = [point_north, point_west]    
+
+point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
+assert zone == refzone
+
+
+#Neils areas around interesting points
+neil1_point1 = [degminsec2decimal_degrees(-20,35,34),
+                degminsec2decimal_degrees(116,45,18)]
+neil1_point2 = [degminsec2decimal_degrees(-20,36,15),
+                degminsec2decimal_degrees(116,46,18)]
+neil1_point3 = [degminsec2decimal_degrees(-20,35,9),
+                degminsec2decimal_degrees(116,47,17)]
+neil1_point4 = [degminsec2decimal_degrees(-20,34,26),
+                degminsec2decimal_degrees(116,46,17)]
+
+neil1_polygon, zone = convert_from_latlon_to_utm([neil1_point1,
+                                                  neil1_point2,
+                                                  neil1_point3,
+                                                  neil1_point4])
+assert zone == refzone
+
+
+
+
+neil2_point1 = [degminsec2decimal_degrees(-20,39,36),
+                degminsec2decimal_degrees(116,41,33)]
+neil2_point2 = [degminsec2decimal_degrees(-20,40,10),
+                degminsec2decimal_degrees(116,42,13)]
+neil2_point3 = [degminsec2decimal_degrees(-20,38,39),
+                degminsec2decimal_degrees(116,43,49)]
+neil2_point4 = [degminsec2decimal_degrees(-20,38,5),
+                degminsec2decimal_degrees(116,43,9)]
+
+neil2_polygon, zone = convert_from_latlon_to_utm([neil2_point1,
+                                                  neil2_point2,
+                                                  neil2_point3,
+                                                  neil2_point4])
+assert zone == refzone
+
+
+
+
+
+#Withnell bay
+wb_point1 = [degminsec2decimal_degrees(-20,35,34),
+                degminsec2decimal_degrees(116,45,18)]
+wb_point2 = [degminsec2decimal_degrees(-20,36,15),
+                degminsec2decimal_degrees(116,46,18)]
+wb_point3 = [degminsec2decimal_degrees(-20,35,9),
+                degminsec2decimal_degrees(116,47,17)]
+wb_point4 = [degminsec2decimal_degrees(-20,34,26),
+                degminsec2decimal_degrees(116,46,17)]
+
+wb_polygon, zone = convert_from_latlon_to_utm([wb_point1, wb_point2,
+                                               wb_point3, wb_point4])
+assert zone == refzone
+
+
+
+
+
+#Larger Withnell bay
+lwb_point1 = [degminsec2decimal_degrees(-20,35,59),
+                degminsec2decimal_degrees(116,42,00)]
+lwb_point2 = [degminsec2decimal_degrees(-20,36,50),
+                degminsec2decimal_degrees(116,46,50)]
+lwb_point3 = [degminsec2decimal_degrees(-20,34,00),
+                degminsec2decimal_degrees(116,47,39)]
+lwb_point4 = [degminsec2decimal_degrees(-20,33,00),
+                degminsec2decimal_degrees(116,42,50)]
+
+lwb_polygon, zone = convert_from_latlon_to_utm([lwb_point1, lwb_point2,
+                                                lwb_point3, lwb_point4])
+                                                      
+assert zone == refzone
+
+
+
+
