#! /usr/bin/env python
# This is based on the Bhutan version.
# This is part of a package for upscaling DEM flowdirections and
# flowaccumulations.  It will be hardwired at 100:1 for 
# development.  It reads the bil file for the buffered mask (at low res)
# for the basin, determines the coordinates for basin cells, and
# assembles the 3*calex3*scale arrays of flowdirection
# for the cell and its neighborhood.
# They represent 9 composite cells.  We calculate how much flow goes
# from the center cell to each of it's neighbors
# The array function supports one-dimensional arrays, so we keep track of
# rows and columns.
# This version reads 2 big bil file: flowdir and flowacc 
# This was originally called readmaster.py.  After ninecells.py was folded
# in, the name "flowscale" seemed more appropriate.
# Always remember to open binary files in binary mode.  Unix forgives, but not XP
# flowscale.py v. 0.84 handles sinks
# Hardwired variables:
# other filenames, ratio in maxlook, diagnostics

# Input files:
#  root.bil		# buffered low-res mekong mask
#  flowdir.bil		# flowdirection, high-res, same bounds
#  flowacc.flt		# flowaccumulation, high-res, same bounds

# Output files: 
#  'blockdir.bil'	# This is what we want.
#  'Errorfile'
#  a8			# diagnostic dump of row,col,tots[]
#  'arrows.gen','w'	# diagnostic file of flow arrows
#  'pointaccs.txt','w' 	# diagnostic file of accs of edge cells
# If the second argument is an extra extension, it will be appended to
#    root and blockdir.  This allows runs with different masks or
#    resolutions in the same directory.
import sys
import os
import string
import array

print 'flowscale.py v. 0.85 introduces arguments: scale, extra extension and infile'
if len(sys.argv) > 1:
  scale=int(sys.argv[1])
else:
  scale = 100 # scale from 3 seconds to 5 minutes
  print 'Using default scaling of 100 (scale from 3 seconds to 5 minutes)'
if len(sys.argv) > 2:
  ext=sys.argv[2]
else:
  ext=""
if len(sys.argv) > 3:
  root=sys.argv[3]
else:
  print 'Using default input filename vicmaskplus8'
  root='vicmaskplus8'	# The one-bit mask was more trouble than it was worth.
			# The mask but be in the aggregate scale, and
			# extend at least one aggregate cell beyond the
			# final basin.
root=root+ext
maxlook=int(.58*scale)	# number of cells to check to see if it goes diagonal
			# was .62 in version 0.81
diagrows = [3,4,5]	# rows for verbose diagnostics
maxdiagnostic = max(diagrows)+1
uplist=[32,64,128]
leftlist=[8,16,32]
downlist=[2,4,8]
rightlist=[1,2,128]
dx=[0]*256	#129
dx[1]=1
dx[2]=1
dx[4]=0
dx[8]=-1
dx[16]=-1
dx[32]=-1
dx[64]=0
dx[128]=1
dy=[0]*256	#129
dy[1]=0
dy[2]=1	# in rowspace, south is positive
dy[4]=1
dy[8]=1
dy[16]=0
dy[32]=-1
dy[64]=-1
dy[128]=-1
scalesq=scale*scale
crows=3*scale	#crows, in ROWs in a Chunk
ccols=crows
f=open(root+".hdr",'r')
f.readline()		# Assume the particular format that imagegrid9.3 uses
f.readline()		# BIL
nrows = int(string.split(f.readline())[1])
ncols = int(string.split(f.readline())[1])
f.readline()
f.readline()
f.readline()
rowbytes = int(string.split(f.readline())[1])
f.close()

f=open(root+".blw",'r')
xdim = float(string.split(f.readline())[0])
f.readline()
f.readline()
f.readline()
ulxmap = float(string.split(f.readline())[0])
ulymap = float(string.split(f.readline())[0])
print nrows,ncols,rowbytes,ulxmap,ulymap
f.close()
xmin=ulxmap-(xdim/2.)
ymax=ulymap+(xdim/2.)

e=open('Errorfile','w')		# Error file for strange file problems
m=open(root+".bil",'rb')		# buffered low-res mekong mask
d=open('flowdir.bil','rb')	# flowdirection, high-res, same bounds
a=open('flowacc.flt','rb')	# flowaccumulation, high-res, same bounds
a8=open('a8','w')		# diagnostic dump of row,col,tots[]
g=open('arrows.gen','w')	# diagnostic file of flow arrows
pt=open('pointaccs.txt','w')	# diagnostic file of accs of edge cells
pt.write('ID\tX\tY\tACC\n')
blockdir=array.array('B',[0]*(nrows*ncols))
blockacc=array.array('f',[0.]*(nrows*ncols))
alltots=[[0,0,0,0,0,0,0,0]]*(nrows*ncols)	# list of lists
mask=[]		# mask in blockspace of cells to process
for row in range(nrows):	# in blockspace
  print  row
  if row == maxdiagnostic:
    g.write('end\n')
    pt.close()
    g.close()
  flags=array.array('B')
  flags.fromfile(m,ncols)
  mask.append(flags)		# odd programming, mixing rows and 2D
  for col in range(ncols):	# in blockspace
    if flags[col]:		# extract 3x3 neighborhood
      ulrow = (row-1)*scale	# in finespace
      ulcol = (col-1)*scale
      dirs=array.array('B')
      accs=array.array('f')
#     print row,col,ulrow,ulcol,0,4*((ulrow+0)*ncols*scale + ulcol)
      for r in range(3*scale):
        d.seek((ulrow+r)*ncols*scale + ulcol)
        dirs.fromfile(d,3*scale)
        seeking = 4*((ulrow+r)*ncols*scale + ulcol)
        a.seek(seeking)
        accs.fromfile(a,3*scale)
        if a.tell() != seeking+12*scale:
          e.write('Bizarre file Error, %d not equal to %d\n' %(a.tell(),seeking+12*scale))
#     print accs[30scale],accs[30101],accs[30400],accs[30401],accs[30700],accs[30701] # diag
#     print dirs[30scale],dirs[30101],dirs[30400],dirs[30401],dirs[30700],dirs[30701] # diag
    # do the work  on the 3*scalex3*scale block here
    # This bit of code started as y:/topog/areas/hydrosheds/5min/ninecells.py
#     tots=array.array('f',[0,0,0,0,0,0,0,0]) 
      tots=[0]*8	# total flows to neighboring cells
      			# number clockwise with eastward=tots[0]
# repetitive code blocks for each corner and each edge
# 1. UL corner
      p=scalesq*3+scale	# upper left cell of center (pointer to 1D array)
      rx=0	# relative x in fine cellspace
      ry=0	# relative y in fourth quadrant
      acc=accs[p]	# flow leaving(?) block
      if row in diagrows:
        pt.write('1\t%f\t%f\t%.0f\n' % (xmin+(col*xdim)+(xdim*0.5/scale),ymax-(row*xdim)-(0.5*xdim/scale),acc))
      if acc > 0:
        if row in diagrows:
          g.write('1\n%f %f\n' % (xmin+(col*xdim)+(xdim*0.5/scale),ymax-(row*xdim)-(0.5*xdim/scale)))
        for j in range(maxlook):	# maxlook chances to find diagonal block
            dir=dirs[p]
            rx=rx+dx[dir]
            ry=ry+dy[dir]
            if row in diagrows:
              g.write('%f %f\n' % (xmin+(col*xdim)+(xdim*(rx+0.5)/scale),ymax-(row*xdim)-((ry+0.5)*xdim/scale)))
            done=0
#           print 'ul',row,col,' rx',rx,', ry',ry,', dir',dir,',acc',acc,', p',p
            if (rx >= 0) & (ry >= 0):	# re-enters block from left
              done=1
              break
            if (rx < 0) & (ry < 0):	# enters upperleft block 
              tots[5] += acc	# flows to upper left cell
              done=1
              break
            p=p+dx[dir]+3*scale*dy[dir]
        if row in diagrows:
          g.write('end\n')
        if not done:
          if rx >= 0:		# stays in upper block
            tots[6] += acc
          elif ry >= 0:
            tots[4] += acc	# stays in left block
          else:
            print 'UL trouble',row,col,done,p,rx,ry
#       print 'UL',done,p,rx,ry,tots
# 2. UR corner
      p=scalesq*3+2*scale-1	# upper right cell of center 
      rx=0	# relative x
      ry=0	# relative y in fourth quadrant
      acc=accs[p]	# flow leaving(?) block
      if row in diagrows:
        pt.write('2\t%f\t%f\t%.0f\n' % (xmin+((col+1)*xdim)-(xdim*0.5/scale),ymax-(row*xdim)-(0.5*xdim/scale),acc))
      if acc > 0:
        for j in range(maxlook):	# maxlook chances to find diagonal block
            dir=dirs[p]
            rx=rx+dx[dir]
            ry=ry+dy[dir]
            done=0
#           print 'ur',row,col,' rx',rx,', ry',ry,', dir',dir,',acc',acc,', p',p
            if (rx <= 0) & (ry >= 0):	# re-enters block from right
              done=1
              break
#           if (rx <= 0) & (ry == 0): 	# re-enters block from top
#             done=1
#             break
            if (rx > 0) & (ry < 0):	# enters upperright block 
              tots[7] += acc	# flows to upper left cell
              done=1
              break
            p=p+dx[dir]+3*scale*dy[dir]
        if not done:
          if rx <= 0:		# stays in upper block
            tots[6] += acc
          elif ry >= 0:
            tots[0] += acc	# stays in right block
          else:
            print 'UR trouble',row,col,done,p,rx,ry
#       print 'UR',done,p,rx,ry,tots
# 3. LR corner
      p=scalesq*6-scale-1	# lower right cell of center 
      rx=0	# relative x
      ry=0	# relative y in fourth quadrant
      acc=accs[p]	# flow leaving(?) block
      if row in diagrows:
        pt.write('3\t%f\t%f\t%.0f\n' % (xmin+((col+1)*xdim)-(xdim*0.5/scale),ymax-((row+1)*xdim)+(0.5*xdim/scale),acc))
      if acc > 0:
        for j in range(maxlook):	# maxlook chances to find diagonal block
            dir=dirs[p]
            rx=rx+dx[dir]
            ry=ry+dy[dir]
            done=0
#           print 'lr',row,col,' rx',rx,', ry',ry,', dir',dir,',acc',acc,', p',p
            if (rx <= 0) & (ry <= 0):	# re-enters block from right
              done=1
              break
            if (rx > 0) & (ry > 0):	# enters lowerright block 
              tots[1] += acc	# flows to upper left cell
              done=1
              break
            p=p+dx[dir]+3*scale*dy[dir]
        if not done:
          if rx <= 0:		# stays in lower block
            tots[2] += acc
          elif ry <= 0:
            tots[0] += acc	# stays in right block
          else:
            print 'LR trouble',row,col,done,p,rx,ry
#       print 'LR',done,p,rx,ry,tots
# 4. LL corner
      p=scalesq*6-2*scale	# lower left cell of center 
      rx=0	# relative x
      ry=0	# relative y in fourth quadrant
      acc=accs[p]	# flow leaving(?) block
      if row in diagrows:
        pt.write('4\t%f\t%f\t%.0f\n' % (xmin+(col*xdim)+(xdim*0.5/scale),ymax-((row+1)*xdim)+(0.5*xdim/scale),acc))
      if acc > 0:
        for j in range(maxlook):	# maxlook chances to find diagonal block
            dir=dirs[p]
            rx=rx+dx[dir]
            ry=ry+dy[dir]
            done=0
#           print 'll',row,col,' rx',rx,', ry',ry,', dir',dir,',acc',acc,', p',p
            if (rx >= 0) & (ry <= 0):	# re-enters block from left
              done=1
              break
            if (rx < 0) & (ry > 0):	# enters lowerleft block 
              tots[3] += acc	# flows to upper left cell
              done=1
              break
            p=p+dx[dir]+3*scale*dy[dir]
        if not done:
          if rx >= 0:		# stays in lower block
            tots[2] += acc
          elif ry <= 0:
            tots[4] += acc	# stays in left block
          else:
            print 'LL trouble',row,col,done,p,rx,ry
#       print 'LL',done,p,rx,ry,tots
# 5. TOP edge
      for x0 in range(1,scale-1):
        p=scalesq*3+scale+x0
        rx=x0	# for #5 we are working relative to upper left
        ry=0
        acc=accs[p]
        if row in diagrows:
          pt.write('5\t%f\t%f\t%.0f\n' % (xmin+(col*xdim)+(xdim*(x0+0.5)/scale),ymax-(row*xdim)-(0.5*xdim/scale),acc))
        if acc > 0:
          if row in diagrows:
            g.write('5\n%f %f\n' % (xmin+(col*xdim)+(xdim*(x0+0.5)/scale),ymax-(row*xdim)-(0.5*xdim/scale)))
          for j in range(maxlook):     # maxlook chances to find diagonal block
            dir=dirs[p]
            rx=rx+dx[dir]
            ry=ry+dy[dir]
            if row in diagrows:
              g.write('%f %f\n' % (xmin+(col*xdim)+(xdim*(rx+0.5)/scale),ymax-(row*xdim)-((ry+0.5)*xdim/scale)))
            done=0
#           print 'top',row,col,'x0',x0,' rx',rx,', ry',ry,', dir',dir,',acc',acc,', p',p
            if rx < 0:
              if ry == 0:	# rare case of diagonal jump
                tots[4] += acc
              else:
                tots[5] += acc	# enters upperleft block
              done=1
              break
            if rx == scale:
              if ry == 0:       # rare case of diagonal jump
                tots[0] += acc
              else:
                tots[7] += acc  # enters upperright block
              done=1
              break
            if ry >= 0:		# re-enter block
              done=1
              break
            p=p+dx[dir]+3*scale*dy[dir]
          if not done:		# stays in top block
            tots[6] += acc
#         print 'TOP',x0,done,p,rx,ry,tots
          if row in diagrows:
            g.write('end\n')

# 6. BOT edge
      for x0 in range(1,scale-1):
        p=scalesq*6-2*scale+x0
        rx=x0	# for #5 we are working relative to lower left
        ry=0
        acc=accs[p]
        if row in diagrows:
          pt.write('6\t%f\t%f\t%.0f\n' % (xmin+(col*xdim)+(xdim*(x0+0.5)/scale),ymax-((row+1)*xdim)+(0.5*xdim/scale),acc))
        if acc > 0:
          if row in diagrows:
            g.write('6\n%f %f\n' % (xmin+(col*xdim)+(xdim*(x0+0.5)/scale),ymax-((row+1)*xdim)+(0.5*xdim/scale)))
          for j in range(maxlook):     # maxlook chances to find diagonal block
            dir=dirs[p]
            rx=rx+dx[dir]
            ry=ry+dy[dir]
            if row in diagrows:
              g.write('%f %f\n' % (xmin+(col*xdim)+(xdim*(rx+0.5)/scale),ymax-((row+1)*xdim)+((0.5-ry)*xdim/scale)))
            done=0
#           print 'bot',row,col,'x0',x0,' rx',rx,', ry',ry,', dir',dir,',acc',acc,', p',p
            if rx < 0:
              if ry == 0:	# rare case of diagonal jump
                tots[4] += acc
              else:
                tots[3] += acc	# enters lowerleft block
              done=1
              break
            if rx == scale:
              if ry == 0:       # rare case of diagonal jump
                tots[0] += acc
              else:
                tots[1] += acc  # enters lowerright block
              done=1
              break
            if ry <= 0:		# re-enter block
              done=1
              break
            p=p+dx[dir]+3*scale*dy[dir]
          if not done:		# stays in bottom block
            tots[2] += acc
#         print 'BOT',x0,done,p,rx,ry,tots
          if row in diagrows:
            g.write('end\n')

# 7. LEFT edge
      for y0 in range(1,scale-1):
        p=scalesq*3+(3*y0+1)*scale
        ry=y0	# for #7 we are working relative to upper left
        rx=0
        acc=accs[p]
        if row in diagrows:
          pt.write('7\t%f\t%f\t%.0f\n' % (xmin+(col*xdim)+(xdim*0.5/scale),ymax-(row*xdim)-((y0+0.5)*xdim/scale),acc))
        if acc > 0:
          if row in diagrows:
            g.write('7\n%f %f\n' % (xmin+(col*xdim)+(xdim*0.5/scale),ymax-(row*xdim)-((y0+0.5)*xdim/scale)))
          for j in range(maxlook):     # maxlook chances to find diagonal block
            dir=dirs[p]
            rx=rx+dx[dir]
            ry=ry+dy[dir]
            if row in diagrows:
              g.write('%f %f\n' % (xmin+(col*xdim)+(xdim*(rx+0.5)/scale),ymax-(row*xdim)-((ry+0.5)*xdim/scale)))
            done=0
#           print 'left',row,col,'y0',y0,' rx',rx,', ry',ry,', dir',dir,',acc',acc,', p',p
            if ry < 0:
              if rx == 0:	# rare case of diagonal jump
                tots[6] += acc
              else:
                tots[5] += acc	# enters upperleft block
              done=1
              break
            if ry == scale:
              if rx == 0:       # rare case of diagonal jump
                tots[2] += acc
              else:
                tots[3] += acc  # enters lowerleft block
              done=1
              break
            if rx >= 0:		# re-enter block
              done=1
              break
            p=p+dx[dir]+3*scale*dy[dir]
          if not done:		# stays in left block
            tots[4] += acc
#         print 'LEFT',y0,done,p,rx,ry,tots
          if row in diagrows:
            g.write('end\n')

# 8. RIGHT edge
      for y0 in range(1,scale-1):
        p=scalesq*3+(3*y0+2)*scale - 1
        ry=y0	# for #8 we are working relative to upper right
        rx=0
        acc=accs[p]
        if row in diagrows:
          pt.write('8\t%f\t%f\t%.0f\n' % (xmin+((col+1)*xdim)-(xdim*0.5/scale),ymax-(row*xdim)-((y0+0.5)*xdim/scale),acc))
        if acc > 0:
          if row in diagrows:
            g.write('8\n%f %f\n' % (xmin+((col+1)*xdim)-(xdim*0.5/scale),ymax-(row*xdim)-((ry+0.5)*xdim/scale)))
          for j in range(maxlook):     # maxlook chances to find diagonal block
            dir=dirs[p]
            rx=rx+dx[dir]
            ry=ry+dy[dir]
            if row in diagrows:
              g.write('%f %f\n' % (xmin+((col+1)*xdim)+(xdim*(rx-0.5)/scale),ymax-(row*xdim)-((ry+0.5)*xdim/scale)))
            done=0
#           print 'right',row,col,'y0',y0,' rx',rx,', ry',ry,', dir',dir,',acc',acc,', p',p
            if ry < 0:
              if rx == 0:	# rare case of diagonal jump
                tots[6] += acc
              else:
                tots[7] += acc	# enters upperright block
              done=1
              break
            if ry == scale:
              if rx == 0:       # rare case of diagonal jump
                tots[2] += acc
              else:
                tots[1] += acc  # enters lowerright block
              done=1
              break
            if rx <= 0:		# re-enter block
              done=1
              break
            p=p+dx[dir]+3*scale*dy[dir]
          if not done:		# stays in right block
            tots[0] += acc
#         print 'RIGHT',y0,done,p,rx,ry,tots
          if row in diagrows:
            g.write('end\n')
      a8.write('%3d %3d  %d %d %d %d %d %d %d %d\n' % (row,col,tots[0],tots[1],tots[2],tots[3],tots[4],tots[5],tots[6],tots[7]))
      cp = (row*ncols)+col	# index to array of coarse cells
#     print row,col,cp,tots
      alltots[cp]=tots
m.close()
d.close()
a.close()
g.close()
maxflows=[0]*nrows*ncols	# info for correcting crisscrosses
for row in range(1,nrows-1):    # in blockspace
  for col in range(1,ncols-1):  # do not start at the edge, lest you flow off
    cp = (row*ncols)+col      # index to array of coarse cells
    if mask[row][col]:
      maxflow=0
      blockdir[cp]=0
# When does a block become a sink?
#	a) Net flow is inward.  (But a block could contain sinks and
#	   pass most of the flow through.
#	b) Total of outflow is less than half the cell.
#	   Does this supercede a) ?  Maybe.
      totnet = 0.
      totout = 0.
      net = alltots[cp][0]-alltots[cp+1][4]        # net flow right
      totnet+=net
      if net > 0.:
        totout+=net
      if net > maxflow:
        maxflow=net
        blockdir[cp] = 1
      net = alltots[cp][1]-alltots[cp+1+ncols][5] # net flow lower right
      totnet+=net
      if net > 0.:
        totout+=net
      if net > maxflow:
        maxflow=net
        blockdir[cp] = 2
      net = alltots[cp][2]-alltots[cp+ncols][6]   # net flow lower cell
      totnet+=net
      if net > 0.:
        totout+=net
      if net > maxflow:
        maxflow=net
        blockdir[cp] = 4
      net = alltots[cp][3]-alltots[cp+ncols-1][7] # net flow lower left
      totnet+=net
      if net > 0.:
        totout+=net
      if net > maxflow:
        maxflow=net
        blockdir[cp] = 8
      net = alltots[cp][4]-alltots[cp-1][0]       # net flow left
      totnet+=net
      if net > 0.:
        totout+=net
      if net > maxflow:
        maxflow=net
        blockdir[cp] = 16
      net = alltots[cp][5]-alltots[cp-ncols-1][1] # net flow upper left
      totnet+=net
      if net > 0.:
        totout+=net
      if net > maxflow:
        maxflow=net
        blockdir[cp] = 32
      net = alltots[cp][6]-alltots[cp-ncols][2]   # net flow upper cell
      totnet+=net
      if net > 0.:
        totout+=net
      if net > maxflow:
        maxflow=net
        blockdir[cp] = 64
      net = alltots[cp][7]-alltots[(cp-ncols)+1][3]  # net flow upper right
      totnet+=net
      if net > 0.:
        totout+=net
      if net > maxflow:
        maxflow=net
        blockdir[cp] = 128
      maxflows[cp]=maxflow
#     if totnet < 0. or totout < scalesq/2:
#       print 'block ',row,col,totnet,totout # do not act yet
      if totout < scalesq/2:
        if totnet < 0. :
          print 'both ',row,col,totnet,totout # do not act yet
        else:
          print 'totout',row,col,totnet,totout # do not act yet
#     These seem to catch sink bottoms.
        blockdir[cp] = 255
#     elif totnet < 0. :
#       print 'totnet ',row,col,totnet,totout # not interesting?


#for row in range(1,nrows-1):    # in blockspace
#  for col in range(1,ncols-1):  # do not start at the edge, lest you flow off
#    cp = (row*ncols)+col      # index to array of coarse cells
#    if mask[row][col]:
#      if blockdir[cp] in [2,8,32,128]:	#diagonal
#        if blockdir[cp]==2:
#          if blockdir[cp+1]==8:
#            if maxflows[cp]<=maxflows[cp+1]:
#              blockdir[cp]=4	# make current cell 4-pointer
#            else:
#              blockdir[cp+1]=4	# make other cell 4-pointer
#          elif blockdir[cp+ncols]==128:	# crisscross the other way
#            if maxflows[cp]<=maxflows[cp+ncols]:
#              blockdir[cp]=1    # make current cell 4-pointer
#            else:
#              blockdir[cp+ncols]=1  # make other cell 4-pointer
#        if blockdir[cp]==8:
#          if blockdir[cp-1]==1:
#            if maxflows[cp]<=maxflows[cp-1]:
#              blockdir[cp]=4	# make current cell 4-pointer
#            else:
#              blockdir[cp-1]=4	# make other cell 4-pointer
#          elif blockdir[cp+ncols]==32:	# crisscross the other way
#            if maxflows[cp]<=maxflows[cp+ncols]:
#              blockdir[cp]=16    # make current cell 4-pointer
#            else:
#              blockdir[cp+ncols]=16  # make other cell 4-pointer
#        if blockdir[cp]==128:
#          if blockdir[cp-ncols]==2:
#            if maxflows[cp]<=maxflows[cp-ncols]:
#              blockdir[cp]=1	# make current cell 4-pointer
#            else:
#              blockdir[cp-ncols]=1	# make other cell 4-pointer
#          elif blockdir[cp+1]==32:	# crisscross the other way
#            if maxflows[cp]<=maxflows[cp+1]:
#              blockdir[cp]=64    # make current cell 4-pointer
#            else:
#              blockdir[cp+1]=64  # make other cell 4-pointer

# from flowscale_mask.py:
for row in range(0,nrows):    # in blockspace
  for col in range(0,ncols):  # do not start at the edge, lest you flow off
    cp = (row*ncols)+col      # index to array of coarse cells
    if basin[row][col]:
      if blockdir[cp] in [2,8,32,128]:  #diagonal
        if blockdir[cp]==2:
          if blockdir[cp+1]==8:
            e.write("crisscross at %d %d, 2 and 8\n" % (row,col))
            if maxflows[cp]<=maxflows[cp+1]:
              blockdir[cp]=4    # make current cell 4-pointer
            else:
              blockdir[cp+1]=4  # make other cell 4-pointer
          elif blockdir[cp+ncols]==128: # crisscross the other way
            e.write("crisscross at %d %d, 2 and 128\n" % (row,col))
            if maxflows[cp]<=maxflows[cp+ncols]:
              blockdir[cp]=1    # make current cell 4-pointer
            else:
              blockdir[cp+ncols]=1  # make other cell 4-pointer
        if blockdir[cp]==8:
          if blockdir[cp-1]==2:
            e.write("Cannot get here at %d %d, 8 and 2\n" % (row,col))
            if maxflows[cp]<=maxflows[cp-1]:
              blockdir[cp]=4    # make current cell 4-pointer
            else:
              blockdir[cp-1]=4  # make other cell 4-pointer
          elif blockdir[cp+ncols]==32:  # crisscross the other way
            e.write("crisscross at %d %d, 8 and 32\n" % (row,col))
            if maxflows[cp]<=maxflows[cp+ncols]:
              blockdir[cp]=16    # make current cell 4-pointer
            else:
              blockdir[cp+ncols]=16  # make other cell 4-pointer
        if blockdir[cp]==32:
          if blockdir[cp-1]==128:
            e.write("Cannot get here at %d %d, 32 and 128\n" % (row,col))
            if maxflows[cp]<=maxflows[cp-1]:
              blockdir[cp]=64    # make current cell 4-pointer
            else:
              blockdir[cp-1]=64  # make other cell 4-pointer
          elif blockdir[cp-ncols]==8:  # crisscross the other way
            e.write("Cannot get here at %d %d, 32 and 8\n" % (row,col))
            if maxflows[cp]<=maxflows[cp-ncols]:
              blockdir[cp]=16    # make current cell 4-pointer
            else:
              blockdir[cp-ncols]=16  # make other cell 4-pointer
        if blockdir[cp]==128:
          if blockdir[cp-ncols]==2:
            e.write("Cannot get here at %d %d, 128 and 2\n" % (row,col))
            if maxflows[cp]<=maxflows[cp-ncols]:
              blockdir[cp]=1    # make current cell 4-pointer
            else:
              blockdir[cp-ncols]=1      # make other cell 4-pointer
          elif blockdir[cp+1]==32:      # crisscross the other way
            e.write("crisscross at %d %d, 128 and 32\n" % (row,col))
            if maxflows[cp]<=maxflows[cp+1]:
              blockdir[cp]=64    # make current cell 4-pointer
            else:
              blockdir[cp+1]=64  # make other cell 4-pointer

f=open('blockdir'+ext+'.bil','wb')
blockdir.tofile(f)
f.close()