source: anuga_work/development/extras/bath2most.py @ 4605

#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()