' **********************************************************
' Gaps_v2.AVE
' Finds gaps and makes new shapefile with them
'
' Modified by Tim Thomas on 12/6/01.
'
' This is from the following:
' Find Gaps and Overlaps
' Author: David M. Theobald, davet@nrel.colostate.edu
' Date: August 1999
' ************************************************************

theView = av.GetActiveDoc
theActiveTheme = theView.GetActiveThemes.Get(0)
theFTab = theActiveTheme.GetFTab
numRecs = theFTab.GetNumRecords
theshapeFld = theFTab.FindField ( "Shape" )
fieldList = theFTab.GetFields.Clone
idField2 = MsgBox.ChoiceAsString(fieldList,"Select ID field",
            "Choose ID field")

' ************************************************************
' Make new shapefile
' ************************************************************

FTabFN = theFTab.GetSrcName.GetFileName
baseString = FTabFN.GetBaseName
extString = FTabFN.GetExtension
newString = FTabFN.AsString.Substitute(baseString,"")
new2String = baseString.Substitute(extString,"shp")
newFTabFN = FileDialog.Put(
          (newString + "gaps " + new2String).AsFileName,
          "*.shp","Gaps shapefile name")
newFTab = FTab.MakeNew(newFTabFN,Polygon) 
newshapeField = newFTab.FindField("Shape")
idField = Field.Make ("GapsID", #FIELD_LONG, 8, 0)
newFTab.AddFields ( { idField } )

' **********************************************************
' Make a comma-delimited file to hold neighbor codes
' **********************************************************

' Shows message in message bar
av.ShowMsg("Creating neighbor file")

' Creates file
new2String = baseString.Substitute(extString,"csv")
neighborFN = FileDialog.Put(
          (newString + "neighbor " + new2String).AsFileName,
          "*.csv","Data file name")
neighborFile = LineFile.Make(neighborFN,#FILE_PERM_WRITE) 


beginNum = MsgBox.Input("Start number for Gaps","Gap number","10000000").AsNumber
kt = 0


' ************************************************************
' FIND GAPS 
' ************************************************************

' Find gaps by going one shape at a time, and finding all shapes within a 
'   rectangle four times the size of the bounding rectangle.  Union them.
'   Then take the difference with the bounding rectangle.

' Shows message in message bar
av.ShowMsg("Step 1: Finding gaps")
av.ShowStopButton

allGaps = theFTab.ReturnValue ( theshapeFld, 0 )

for each rec in theFTab
  thePoly = theFTab.ReturnValue ( theshapeFld, rec )
  theRect = thePoly.ReturnExtent
  x1 = theRect.GetLeft - (theRect.GetWidth * 0.5)
  x2 = theRect.GetRight + (theRect.GetWidth * 0.5)
  y1 = theRect.GetBottom - (theRect.GetHeight * 0.5)
  y2 = theRect.GetTop + (theRect.GetHeight * 0.5)
  bigRect = Rect.MakeXY ( x1, y1, x2, y2 )
  theSel = theFTab.SelectByRect(bigRect,#VTAB_SELTYPE_NEW)
  theUnion = theFTab.ReturnValue ( theshapeFld, 0 )
  for each s in theFTab.GetSelection
    thePoly = theFTab.ReturnValue ( theshapeFld, s )
    theUnion = thePoly.ReturnUnion( theUnion )
  end
  theGap = bigRect.ReturnDifference( theUnion )
  if (theGap.ReturnArea > 0) then
    ' If there is a gap (area > 0), then merge it with other gaps
    allGaps = theGap.ReturnUnion( allGaps )
  end    
  av.SetStatus(100*(kt+1)/numRecs)
  kt = kt + 1
end

gapList = allGaps.Explode
areaList = List.Make
last = gapList.Count - 1
largeArea = 0

if (last > 0) then

  av.ShowMsg("Step 2: Making gaps shapefile")

  for each i in 0..last
    theArea = gapList.Get(i).ReturnArea
    areaList.Add(theArea)
    if (theArea > largeArea) then
      bigID = i
      largeArea = theArea
    end
  end
  
  kt = 0  
  for each i in 0..last
    if (i <> bigID) then
      rec = newFTab.AddRecord
      newFTab.SetValue ( newshapeField, rec, gapList.Get(i) )
      kt = kt + 1
      newFTab.SetValue ( idField, rec, (beginNum - 1 + kt) )
    end
  end
    
  ' ****************************************************************************
  ' The above finishes the work of finding the gaps
  ' The below tries to identify gaps that have more than 2 neighbors.
  '   These are the ones that need more work.  Someone needs to manually extend
  '   the vertices, so that gaps only have 2 neighbors.
  
  av.ShowMsg("Step 3: Checking for neighbors of gaps")
  dwrite = 0
  for each rec in newFTab
    thePoly = newFTab.ReturnValue ( theshapeFld, rec )
    theSel = theFTab.SelectByPolygon(thePoly,#VTAB_SELTYPE_NEW)
    if (theFTab.GetSelection.Count > 2) then
    
      if (dwrite = 0) then
        dwrite = 1
        neighborFile.WriteElt("Data only written when gap has more than 2 neighbors")
        neighborFile.WriteElt("GapID  ShapeID")
      end
      
      id = newFTab.ReturnValueString( idfield, rec )
      
      for each s in theFTab.GetSelection
        ' Writing a line of output to the file
        id2 = theFTab.ReturnValueString( idfield2, s )
        theString = id + ", " + id2
        neighborFile.WriteElt(theString)
      end
      
    end
  end
  
  ' This processes a new shapefile with the gaps combined with the original polys
  '   Only enters this if we have 2 or less neighbors per gap
  if (dwrite = 0) then
  
    ' Make new shapefile
    new2String = baseString.Substitute(extString,"shp")
    new2FTabFN = FileDialog.Put(
          (newString + "no gaps " + new2String).AsFileName,
          "*.shp","Gaps shapefile name")
    new2FTab = FTab.MakeNew(new2FTabFN,Polygon) 
    new2shapeField = new2FTab.FindField("Shape")

    ' Clone fields for newFTab
    inFields = theFTab.GetFields
    newFields = List.Make
    for each f in inFields
      if (f.GetName <> "shape") then
        newFields.Add(f.Clone)
      end
    end
    new2FTab.AddFields(newFields)
    
  
    joinIDList = List.Make
    joinShapeList = List.Make
  
    for each rec in newFTab
      thePoly = newFTab.ReturnValue ( newshapeField, rec )
      theSel = theFTab.SelectByPolygon(thePoly,#VTAB_SELTYPE_NEW)
      kt = 1
      for each s in theFTab.GetSelection
        if (kt = 1) then
          onePoly = theFTab.ReturnValue ( theshapeFld, s )
          oneID = theFTab.ReturnValue( idfield2, s )
          kt = kt + 1
        elseif (kt = 2) then
          twoPoly = theFTab.ReturnValue ( theshapeFld, s )
          twoID = theFTab.ReturnValue( idfield2, s )
        end
      end
        
      ' Pick the polygon to join the gap to (choose the smallest of the two)
      if (theFTab.GetSelection.Count = 1) then
        joinID = oneID
        joinPoly = onePoly
      else
        if (onePoly.ReturnArea < twoPoly.ReturnArea) then
          joinID = oneID
          joinPoly = onePoly
        else
          joinID = twoID
          joinPoly = twoPoly
        end
      end
        
      ' Check to see if the joinPoly has already been joined to a different gap
      '   If so, combine the new gap with the already joined poly
      dfind = 0
      if (joinIDList.Count > 0) then
        ' for each i in joinIDList
        for each i in 0..(joinIDList.Count-1)
          if (joinID = joinIDList.Get(i)) then
            unionPoly = thePoly.ReturnUnion( joinShapeList.Get(i) )
            joinShapeList.Set (i, unionPoly)
            dfind = 1
          end
        end
      end
        
      ' If not already in the joined poly list, add to it
      if (dfind = 0) then
        unionPoly = thePoly.ReturnUnion( joinPoly )
        joinShapeList.Add(unionPoly)
        joinIDList.Add(joinID)
      end
        
    end
    theActiveTheme.ClearSelection
    
    ' After going through all the gaps and joining, we now have to output the 
    '   new dataset
    av.ShowMsg("Step 4: Outputting new file with gaps merged")
    av.ShowStopButton
    kt = 0
    for each rec in theFTab
      thePoly = theFTab.ReturnValue ( theshapeFld, rec )
      theID = theFTab.ReturnValue( idfield2, rec )
      ' Checking to see if shape was joined with a gap
      ' for each i in joinIDList
      for each i in 0..(joinIDList.Count-1)
        if (theID = joinIDList.Get(i)) then
          thePoly = joinShapeList.Get(i)
        end
      end
      newrec = new2FTab.AddRecord
      new2FTab.SetValue(new2shapeField,newrec,thePoly)
      for each fld in theFTab.GetFields
        fldName = fld.GetName
        if (fldName <> "shape") then
          newfld = new2FTab.FindField(fldName)
          val = theFTab.ReturnValue(fld,rec)
          new2FTab.SetValue(newfld,newrec,val)
        end
      end
      av.SetStatus(100*(kt+1)/numRecs)
      kt = kt + 1
    end

    ' finish up the new theme
    new2FTab.SetEditable ( FALSE )
    new2FTheme = FTheme.Make ( new2FTab )
    theView.AddTheme ( new2FTheme )

    MsgBox.Info("All gaps had 2 neighbors or less","")

  end

else

  MsgBox.Info("No gaps were found","No gaps were found")

end

theActiveTheme.ClearSelection 
av.ClearMsg
av.ClearStatus

' finish up the new theme
newFTab.SetEditable ( FALSE )
newFTheme = FTheme.Make ( newFTab )
theView.AddTheme ( newFTheme )
newFTheme.ClearSelection