Index: /anuga_work/development/extras/__init__.py
===================================================================
--- /anuga_work/development/extras/__init__.py	(revision 4540)
+++ /anuga_work/development/extras/__init__.py	(revision 4540)
@@ -0,0 +1,8 @@
+"""Make directory available as a Python package
+"""
+
+pass
+
+
+
+
Index: /anuga_work/development/extras/bath2most.py
===================================================================
--- /anuga_work/development/extras/bath2most.py	(revision 4540)
+++ /anuga_work/development/extras/bath2most.py	(revision 4540)
@@ -0,0 +1,321 @@
+#y_min = -35.07 #100m Woolongong 
+#y_max = -34.+(-10./.6)*0.01
+#x_min = 150.+(+42./.6)*0.01
+#x_max = 151.6
+
+y_min = -35 #100m Botany Bay
+y_max = -33.8
+x_min = 150.+(+50./.6)*0.01
+x_max = 152.4
+
+#y_min = -34.15233 #10m Botany Bay
+#y_max = -33.8950794787
+#x_min = 151.031182
+#x_max = 151.364182675
+
+
+
+#x_min=150+.88# Port Kembla 10m
+#x_max=150+.575/.6
+#y_min=-34.484
+#y_max=-34-0.20/.6
+
+
+#x_min=150+.525/.6 # Port Kembla 50m
+#x_max=150+.575/.6
+#y_min=-34-0.287/.6
+#y_max=-34-0.23/.6
+
+#x_min=150+.525/.6 # Port Kembla
+#x_max=150+.555/.6
+#y_min=-34-0.287/.6
+#y_max=-34-0.27/.6
+
+#x_min=150+.538/.6 # Flagstaff Point
+#x_max=150+.56/.6
+#y_min=-34-0.265/.6
+#y_max=-34-0.245/.6
+
+
+#x_min=150+.538/.6 # Woolongong Beach
+#x_max=150+.555/.6
+#y_min=-34-.27/.6
+#y_max=-34-.26/.6
+
+
+thin_factor = 2
+add_tide = True
+tide_height = 1.1
+clean_bad_points = True
+no_data_value=-1000000
+#if clean_bad_points is False, bad_points will just be 
+#set to no_data_value
+
+crop =True
+#crop =False
+invert_bathymetry = True
+
+input_file = 'dem_geog.asc'
+#output_file = 'wool_200m.txt'
+output_file = 'Botany_Bay_200m.txt'
+
+#input_file = 'wollg_geo_r.asc'
+#output_file = 'Port_Kembla_10m.txt'
+
+#input_file = 'bb_geog.asc'
+#output_file = 'Botany_Bay_10m.txt'
+
+#################################################
+#END OF INPUT
+
+
+
+
+"""
+Takes a bathymetry file of format:
+""
+ncols         1333
+nrows         1030
+xllcorner     151.031182
+yllcorner     34.15233
+cellsize      0.00025000050663948
+NODATA_value  -9999
+56.3 55.36667 54.1 53.03333 52.2 ... (row n)
+54.18333 53.11667 52.26667 51.6 51.33334 ... (row n-1)
+...
+54.18333 53.11667 52.26667 51.6 51.33334 ... (row 1)
+""
+
+and makes it format:
+""
+nx ny (ncols=longitude_range,nrows=latitude_range)
+151.031182(longitude 1)
+151.031182+cellsize(longitude 2)
+151.031182+2*cellsize(longitude 3)
+...
+151.031182+(n-1)*cellsize(longitude n)
+34.1523+cellsize(latitude 1)
+34.1523+2*cellsize(latitude 2)
+34.1523+3*cellsize(latitude 3)
+...
+34.1523+(n-1)*cellsize(latitude n)
+56.3 55.36667 54.1 53.03333 52.2 (row 1)
+54.18333 53.11667 52.26667 51.6 51.33334 (row 2)
+...
+54.18333 53.11667 52.26667 51.6 51.33334 (n)
+
+""
+If llcorner is negative for longitude then it will 
+then it will automatically change the x increment to -cellsize
+
+Options will include cropping and thinning and finding values
+for NANs.
+
+NB it assumes that the system has the memory to cope with any
+give row or column, or the cropped bathymetry, but not nessesarily
+the whole bathymetry.
+"""
+
+thin_factor = int(thin_factor)
+#from Numeric import allclose
+
+in_file = open(input_file,'r')
+
+nx_in = in_file.readline()
+nx_in = nx_in.split()[1]
+nx_in= int(nx_in)
+
+ny_in = in_file.readline()
+ny_in = ny_in.split()[1]
+ny_in = int(ny_in)
+
+x_min_in = in_file.readline()
+x_min_in = x_min_in.split()[1]
+x_min_in = float(x_min_in)
+y_min_in = in_file.readline()
+y_min_in = y_min_in.split()[1]
+y_min_in = float(y_min_in)
+
+cellsize = in_file.readline()
+cellsize = cellsize.split()[1]
+cellsize = float(cellsize)
+
+nodata = in_file.readline()
+nodata = nodata.split()[1]
+nodata = float(nodata)
+
+x_inc_in = cellsize
+y_inc_in = cellsize
+
+x_max_in = x_min_in+(nx_in-1)*x_inc_in
+y_max_in = y_min_in+(ny_in-1)*y_inc_in
+print 'minimum possible x =',x_min_in
+print 'maximum possible x =',x_max_in
+print 'minimum possible y =',y_min_in
+print 'maximum possible y =',y_max_in
+
+x_inc = x_inc_in*thin_factor
+y_inc = y_inc_in*thin_factor
+
+if crop:
+    x_min = x_min_in+x_inc*int((x_min-x_min_in)/x_inc)
+    x_max = x_min_in+x_inc*int((x_max-x_min_in)/x_inc)
+    y_min = y_min_in+y_inc*int((y_min-y_min_in)/y_inc)
+    y_max = y_min_in+y_inc*int((y_max-y_min_in)/y_inc)
+else:
+    x_min = x_min_in
+    x_max = x_max_in
+    y_min = y_min_in
+    y_max = y_max_in
+    x_min = x_min_in+x_inc*int((x_min-x_min_in)/x_inc)
+    x_max = x_min_in+x_inc*int((x_max-x_min_in)/x_inc)
+    y_min = y_min_in+y_inc*int((y_min-y_min_in)/y_inc)
+    y_max = y_min_in+y_inc*int((y_max-y_min_in)/y_inc)
+
+assert not x_min<x_min_in
+assert not y_min<y_min_in
+assert not x_max>x_max_in
+assert not y_max>y_max_in
+assert not x_min>x_max
+assert not y_min>y_max
+
+nx = int(round((x_max-x_min)/x_inc+1))
+ny = int(round((y_max-y_min)/y_inc+1))
+
+x_max = x_min+(nx-1)*x_inc
+y_max = y_min+(ny-1)*y_inc
+
+print 'x_min =',x_min
+print 'x_max =',x_max
+print 'y_min =',y_min
+print 'y_max =',y_max
+
+x_min_index = int(round((x_min-x_min_in)/x_inc_in))
+x_max_index = int(round((x_max-x_min_in)/x_inc_in))
+
+y_min_index = int(round((y_max_in-y_min)/y_inc_in))
+y_max_index = int(round((y_max_in-y_max)/y_inc_in))
+
+assert x_max==x_min+(nx-1)*x_inc
+assert y_max==y_min+(ny-1)*y_inc
+
+assert x_max_in==x_min_in+(nx_in-1)*x_inc_in
+assert y_max_in==y_min_in+(ny_in-1)*y_inc_in
+
+#assert allclose(x_min_index*x_inc_in+x_min_in,x_min)
+#assert allclose(y_min_index*y_inc_in+y_min_in,y_min)
+#assert abs(x_min_index*x_inc_in+x_min_in-x_min)<0.0000001
+#assert abs(y_min_index*y_inc_in+y_min_in-y_min)<0.0000001
+
+#assert allclose(x_max_index*x_inc_in+x_min_in,x_max)
+#assert allclose(y_max_index*y_inc_in+y_min_in,y_max)
+#assert abs(x_max_index*x_inc_in+x_min_in-x_max)<0.0000001
+#assert abs(y_max_index*y_inc_in+y_min_in-y_max)<0.0000001
+
+#assert abs((y_min_index*1.)/(ny_in-1)-(y_min-y_min_in)/(y_max_in-y_min_in))<0.0000001
+
+h1_list=[]
+for i in range(nx):
+    h1_list.append(x_min+i*x_inc)
+
+h2_list=[]
+for j in range(ny):
+    h2_list.append(y_min+j*y_inc)
+
+h2_list.reverse()
+
+print 'reading bathymetry'
+#burn till line = first wanted line
+for i in range(y_max_index):
+    line = in_file.readline()
+
+
+if invert_bathymetry:
+    up=-1.
+
+nodata=nodata*up+tide_height*up*add_tide*-1.
+
+out_depth_array = []
+nodata_dict = {}
+for i in range(y_min_index-y_max_index+1):
+    k=0
+    in_line = in_file.readline()
+    if (i/thin_factor)*thin_factor-i == 0:
+        print 'reading line %i'%(i+y_min_index)
+        print len(in_line.split())
+        out_depth_array.append([])
+        for string in in_line.split()[x_min_index:x_max_index+1]:
+            if (k/thin_factor)*thin_factor-k == 0:
+                out_depth_array[i/thin_factor].append\
+                  (float(string)*up+tide_height*up*add_tide*-1.)
+                if out_depth_array[i/thin_factor][k/thin_factor]==nodata:
+                    nodata_dict[(i/thin_factor,k/thin_factor)]=None
+            k+=1
+
+in_file.close()
+
+
+if clean_bad_points:
+    print 'cleaning %i nodata points'%len(nodata_dict)
+#cleaning NANs
+    bad_points = {}
+    for point in nodata_dict.keys():
+        bad_points[point]=None
+#
+    for j in range(-1,nx+1):
+        point = (-1,j)
+        bad_points[point]=None
+        point = (ny,j)
+        bad_points[point]=None
+#
+    for i in range(-1,ny+1):
+        point = (i,-1)
+        bad_points[point]=None
+        point = (i,nx)
+        bad_points[point]=None
+#
+    while len(nodata_dict)>0:
+        for point in nodata_dict.keys():
+            sum = 0
+            n = 0
+            for i in (-1+point[0],0+point[0],1+point[0]):#for neighboring points
+                for j in (-1+point[1],0+point[1],1+point[1]):
+                    if not bad_points.has_key((i,j)):
+                        sum+=out_depth_array[i][j]
+                        n+=1
+            if n>0:
+                out_depth_array[point[0]][point[1]]=sum/n
+                bye = nodata_dict.pop(point)
+                bye = bad_points.pop(point)
+        print '%i points to go...'%len(nodata_dict)
+
+if not clean_bad_points:
+     for point in nodata_dict.keys():
+         out_depth_array[point[0]][point[1]]=no_data_value
+
+#out_depth_array.reverse()
+
+print 'writing results'
+out_file = open(output_file,'w')
+
+out_file.write('%i %i\n'%(nx,ny))
+
+for line in h1_list:
+    out_file.write('%10.10f'%line)
+    out_file.write('\n')
+
+for line in h2_list:
+    out_file.write('%10.10f'%line)
+    out_file.write('\n')
+
+k=0
+for line in out_depth_array:
+#    line.reverse()
+    print 'writing line %i'%k
+    for depth in line:
+        out_file.write('%10.10f '%depth)
+    out_file.write('\n')
+    k+=1
+
+
+out_file.close()
Index: /anuga_work/development/extras/m2ferret.py
===================================================================
--- /anuga_work/development/extras/m2ferret.py	(revision 4540)
+++ /anuga_work/development/extras/m2ferret.py	(revision 4540)
@@ -0,0 +1,99 @@
+#input_file = 'wool_100m.txt'
+#output_file = 'wool_100m_fer.txt'
+#variable = 'wool_100m'
+#journal_file = 'wool_100m.jnl'
+#journal_file = 'Port_Kembla_10m.jnl'
+#variable = 'Port_Kembla_10m'
+#input_file = 'Port_Kembla_10m.txt'
+#output_file = 'Port_Kembla_10m_fer.txt'
+#input_file = 'Botany_Bay_10m.txt'
+#output_file = 'Botany_Bay_fer_10m.txt'
+#variable = 'Botany_Bay_10m'
+#journal_file = 'Botany_Bay_10m.jnl'
+
+#input_file = 'Botany_Bay_100m.txt'
+#output_file = 'Botany_Bay_fer_100m.txt'
+#variable = 'Botany_Bay_100m'
+#journal_file = 'Botany_Bay_100m.jnl'
+
+#input_file = 'Botany_Bay_200m.txt'
+#output_file = 'Botany_Bay_fer_200m.txt'
+#variable = 'Botany_Bay_200m'
+#journal_file = 'Botany_Bay_200m.jnl'
+
+input_file = 'wool_200m.txt'
+output_file = 'wool_200m_fer.txt'
+variable = 'wool_200m'
+journal_file = 'wool_200m.jnl'
+
+#input_file = 'aus_coast_most_smooth.txt'
+#output_file = 'aus_coast_fer_smooth.txt'
+
+in_file = open(input_file,'r')
+print 'reading header'
+nx_ny_str = in_file.readline()
+nx_str,ny_str = nx_ny_str.split()
+nx = int(nx_str)
+ny = int(ny_str)
+h1_list=[]
+for i in range(nx):
+    h1_list.append(in_file.readline())
+
+x_min = float(h1_list[0])
+x_max = float(h1_list[-1])
+x_inc = (x_max-x_min)*1.0000001/float(nx-1)
+
+h2_list=[]
+for j in range(ny):
+    h2_list.append(in_file.readline())
+
+y_min = float(h2_list[-1])
+y_max = float(h2_list[0])
+y_inc = (y_max-y_min)*1.00000001/float(ny-1)
+
+print 'reading depths'
+in_depth_list = in_file.readlines()
+in_file.close()
+
+out_depth_list = []
+
+for j in range(ny):
+    out_depth_list.append('')
+
+print 'processing depths'
+k=1
+for in_line in in_depth_list:
+    for string in in_line.split():
+        #j = k/nx
+        out_depth_list[(k-1)/nx]=out_depth_list[(k-1)/nx]+('%s    '%string)
+        if k-(k/nx)*nx ==0:
+            out_depth_list[(k-1)/nx]=out_depth_list[(k-1)/nx]+'\n'
+        if k==nx*ny:
+            break
+        k+=1
+
+in_file.close()
+
+out_depth_list.reverse()
+
+depth_list = out_depth_list
+
+print 'writing results'
+out_file = open(output_file,'w')
+
+for line in depth_list:
+    out_file.write(line)
+
+out_file.close()
+
+
+jnl_file = open(journal_file,'w')
+
+jnl_file.write('DEFINE AXIS/X=%10.10fe:%10.10fe:%10.10f/unit=degree xlong\n'%(x_min,x_max,x_inc))
+jnl_file.write('DEFINE AXIS/Y=%10.10fn:%10.10fn:%10.10f/unit=degree ylat\n'%(y_min,y_max,y_inc))
+jnl_file.write('DEFINE GRID/x=xlong/y=ylat Pacific\n')
+jnl_file.write('FILE/VARIABLES=%s/COLUMNS=%i/GRID=Pacific %s\n'%(variable,nx,output_file))
+jnl_file.write('shade/pal=land_sea/lev=(-60,60,1) -1*%s\n'%variable)
+jnl_file.close()
+
+ 
Index: /anuga_work/development/extras/thinbathy.py
===================================================================
--- /anuga_work/development/extras/thinbathy.py	(revision 4540)
+++ /anuga_work/development/extras/thinbathy.py	(revision 4540)
@@ -0,0 +1,49 @@
+"""Thin out MOST NetCDF bathymetry and store as nc file
+
+"""
+
+#input_file = 'e5_bathy_corr5.cdf'
+#output_file = 'SU-AU_e.nc'
+
+#This should have a test
+
+def thinbathy(input_file, output_file):
+
+    lat_thin = lon_thin = 4 
+
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    fid = NetCDFFile(input_file, 'r')
+    print fid.variables
+
+    # for first event
+    latitudes = fid.variables['LAT'][::lat_thin]
+    longitudes = fid.variables['LON'][::lon_thin]
+    elevations = fid.variables['ELEVATION'][::lat_thin,::lon_thin]
+
+    outfile = NetCDFFile(output_file, 'w')
+    outfile.createDimension('LON', fid.dimensions['LON']/lon_thin)
+    outfile.createDimension('LAT', fid.dimensions['LAT']/lat_thin)
+    
+##    # for second event
+##    latitudes = fid.variables['LAT601_690'][::lat_thin]
+##    longitudes = fid.variables['LON645_705'][::lon_thin]
+##    elevations = fid.variables['ELEVATION'][::lat_thin,::lon_thin]
+##
+##    outfile = NetCDFFile(output_file, 'w')
+##    outfile.createDimension('LON', fid.dimensions['LON645_705']/lon_thin + 1)
+##    outfile.createDimension('LAT', fid.dimensions['LAT601_690']/lat_thin + 1)
+
+    import Numeric
+    outfile.createVariable('LON', Numeric.Float64, ('LON',))
+    outfile.createVariable('LAT', Numeric.Float64, ('LAT',))
+    outfile.createVariable('ELEVATION', Numeric.Float32, ('LAT', 'LON'))
+
+    #Assign to NC file
+    #Precision used by MOST for lat/lon is 4 or 5 decimals
+    outfile.variables['LAT'][:] = Numeric.around(latitudes, 5)
+    outfile.variables['LON'][:] = Numeric.around(longitudes, 5)
+    outfile.variables['ELEVATION'][:] = elevations
+
+    fid.close()
+    outfile.close()
