' Written by Timothy S. Thomas.
' Programming began 11/18/99.
' Last update was 12/9/99.


' ***********************************************************************
' This script computes areas of polygons for themes in lat-long, by first
' projecting them into an equal area projection.  It then puts the 
' computed areas into a dbf file, and joins them to the original
' shapefile.
' ***********************************************************************

' Beginning program
theProject = av.GetProject
theView = av.GetActiveDoc
theTheme = theView.GetActiveThemes.Get(0)
theFTab = theTheme.GetFTab


' **********************************************************
' Establish unique polygon id
' **********************************************************

' This allows the user to pick a field as the basic unit for
' evaluation.
fieldList = theFTab.GetFields.Clone
fieldList.Insert("<Make new unique id field>")
idString = MsgBox.ChoiceAsString(fieldList,"Select ID field",
            "Choose ID field")
if (idString = "<Make new unique id field>") then
  theFTab.SetEditable(True)
  isFld = theFTab.FindField("UniqID")
  if (isFld <> nil) then
    idnew = MsgBox.YesNo("UniqID is already a field.  Do you want to use it?",
                         "Use UniqID?",True)
    if (idnew) then
      idField = isFld
    else
      newName = MsgBox.Input("Name of new ID","New name","NewID")
      idField = Field.Make(newName,#FIELD_LONG,8,0)
      theFTab.AddFields({idField})
      for each rec in theFTab
        theFTab.SetValue(idField,rec,rec+1)
      end
    end
  else
    idField = Field.Make("UniqID",#FIELD_LONG,8,0)
    theFTab.AddFields({idField})
    for each rec in theFTab
      theFTab.SetValue(idField,rec,rec+1)
    end
  end
  theFTab.SetEditable(False)
else
  idField = idString
end
idString = idField.AsString

' **********************************************************
' Change projection of new view
' **********************************************************

' Rectangle for projection
theRect = theTheme.ReturnExtent

' Ask which projection
prjList = {"Albers Equal Area Conic","Cylindrical Equal Area",
  "Lambert Azimuthal Equal Area"}
prjName = MsgBox.ChoiceAsString(prjList,"Select an equal area projection",
  "Select a Projection")
 
' Ask units for area
unitList = {"Square kilometers","Square meters","Hectares"} 
unitName = MsgBox.ChoiceAsString(unitList,"Select units for area",
  "Select Units")

' Make projection  
if (prjName = "Albers Equal Area Conic") then
  thePrj = Albers.Make(theRect)
  lowpar = MsgBox.Input("Type the lower standard parallel, or <Enter> for default calculation",
          "Lower Standard Parallel",(5*theRect.GetBottom + theRect.GetTop)/6.AsString)
  thePrj.SetLowerStandardParallel(lowpar.AsNumber)
  uppar = MsgBox.Input("Type the upper standard parallel, or <Enter> for default calculation",
          "Upper Standard Parallel",(theRect.GetBottom + 5*theRect.GetTop)/6.AsString)
  thePrj.SetUpperStandardParallel(uppar.AsNumber)
elseif (prjName = "Cylindrical Equal Area") then
  thePrj = EqualAreaCylindrical.Make(theRect)
elseif (prjName = "Lambert Azimuthal Equal Area") then
  thePrj = EqualAreaAzimuthal.Make(theRect)
end

' Set parameters for projection
merid = MsgBox.Input("Type central meridian, or <Enter> for default of center of theme",
          "Longitude (X)",theRect.ReturnCenter.GetX.AsString)
thePrj.SetCentralMeridian(merid.AsNumber)
lat = MsgBox.Input("Type reference latitude, or <Enter> for default of center of theme",
          "Latitude (Y)",theRect.ReturnCenter.GetY.AsString)
thePrj.SetReferenceLatitude(lat.AsNumber)
  
if (unitName = "Square kilometers") then
  theView.SetUnits(#UNITS_LINEAR_KILOMETERS)
  theView.GetDisplay.SetDistanceUnits(#UNITS_LINEAR_KILOMETERS)
else
  theView.SetUnits(#UNITS_LINEAR_METERS)
  theView.GetDisplay.SetDistanceUnits(#UNITS_LINEAR_METERS)
end

theView.SetProjection(thePrj)

prjFN = "Temp01.shp".AsFileName
prjFTab = theFTab.ExportProjected(prjFN,thePrj,theFTab.GetSelection.Count > 0)
    ' The last term says that if nothing is selected, project all; otherwise, project
    ' the selection.
prjTheme = FTheme.Make(prjFTab)
newView = View.Make
newView.GetWin.Open
newView.AddTheme(prjTheme)
nullPrj = Prj.MakeNull
theView.SetProjection(nullPrj)

' A nice way to get a directory name and use it
FTabFN = theFTab.GetSrcName.GetFileName
baseString = FTabFN.GetBaseName
extString = FTabFN.GetExtension
newString = FTabFN.AsString.Substitute(baseString,"")
new2String = baseString.Substitute(extString,"dbf")
theFN = FileDialog.Put(
          (newString + "Area of " + new2String).AsFileName,
          "*.dbf","Data file name")
dataVTab = VTab.MakeNew(theFN,dBase)

prjidField = prjFTab.FindField(idString)
newidField = prjidField.Clone
newAreaString = MsgBox.Input("8 characters or less, no spaces",
                             "Variable name for area field","NewArea")
newAreaField = Field.Make(newAreaString,#Field_Double,16,2)
dataVTab.AddFields({newidField,newAreaField})

dataVTab.SetEditable(True)
prjFTab.SetEditable(True)
shapeField = prjFTab.FindField("Shape")
numRecs = prjFTab.GetNumRecords

' Shows message in message bar
av.ShowMsg("Calculating area")

for each r in prjFTab
  thePoly = prjFTab.ReturnValue(shapeField,r)
  ' thePolyDense = thePoly.ReturnDensified(5)
  ' thePoly = Polygon.Make(thePolyDense.AsList)
  ' prjFTab.SetValue(shapeField,r,thePoly)
  outrec = dataVTab.AddRecord
  dataVTab.SetValue(newidField,outrec,prjFTab.ReturnValue(prjidField,r))
  if (unitName = "Square kilometers") then
    dataVTab.SetValue(newAreaField,outrec,thePoly.ReturnArea)
  elseif (unitName = "Square meters") then
    dataVTab.SetValue(newAreaField,outrec,thePoly.ReturnArea)
  else  
    dataVTab.SetValue(newAreaField,outrec,thePoly.ReturnArea/10000)
  end
  av.SetStatus(100*(r+1)/numRecs)
end 

' Clears message in message bar
av.ClearMsg

removeList = List.Make
for each d in theProject.GetDocs
  if (d.AsString = newAreaString and d.Is(Table)) then
    removeList.Add(d)
  end
end
for each d in removeList
  theProject.RemoveDoc(d)
end
dataTable = Table.Make(dataVTab)
dataTableWin = dataTable.GetWin
dataTableWin.Open
dataTable.SetName(newAreaString)
theFTable = theTheme.EditTable

' Join DBF file to FTab
theFTab.Join(idField, dataVTab, newidField)

' Turns off editing
theFTab.SetEditable(false)
dataVTab.SetEditable(false)

dataTableWin.Close


newView.DeleteTheme(prjTheme)
av.GetProject.RemoveDoc(newView)
av.PurgeObjects