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