[postgis-users] How to generalize or simplify a Polygon

[oscar.gomez at telvent.abengoa.com]

This is an ArcView-Avenue implementation of the Douglas - Peucker 
algorithm to generalize a polyline or polygon feature's vertices.

Good luck,

Oscar Gomez
IT Consultant
Telvent


theView = av.GetActiveDoc
thePrj = theView.GetProjection
if (thePrj.IsNull) then
  hasPrj = false
else
  hasPrj = true
  thePrj = theView.GetProjection
end

' Get active themes

theThemes = theView.GetActiveThemes

if (theThemes.Count < 1) then
  MsgBox.Error("Please make a polyline or polygon theme active","*** Error 
***")
  exit
end

' Process each active theme

for each t in theThemes

  ' Check that the theme is of the correct type

  if (t.Is( FTHEME).Not) then
    MsgBox.Error (t.GetName++"is not a polyline or polygon theme","*** 
Error ***")
    continue
  end

  theClassName = t.GetFTab.GetShapeClass.GetClassName
  if ((theClassName = "Polygon").Not and (theClassName = "Polyline").Not) 
then
    MsgBox.Error (t.GetName++"is not a polyline or polygon theme","*** 
Error ***")
    continue
  end

  ' Get a tolerance in (distance units) to be used (must be > 0)

  theTol = 0
  theMUnits = theView.GetDisplay.GetUnits
  theDUnits = theView.GetDisplay.GetDistanceUnits
  theUStr = Units.GetFullUnitString(theDUnits).LCase
  while ((theTol <= 0) or (theTol.AsString.IsNumber).Not)
    theTol = MsgBox.Input("Positive tolerance value (in "+theUStr+
         ") for theme"++t.GetName,"Tolerance",theTol.AsString)
    if (theTol = NIL) then exit end

    if (theTol.IsNumber) then
      theTol = theTol.AsNumber
    end
  end
  theTol = Units.Convert(theTol,theDUnits,theMUnits)

  ' Get the name for a new theme, from View.Export script

  def = av.GetProject.MakeFileName("theme", "shp")
  def = FileDialog.Put(def, "*.shp", "New file for simplified "+ 
t.getName)

  if (def = NIL) then return nil end

  ' Create a new shapefile theme, from View.Export script

  tbl = t.GetFtab
  shpfld = (tbl.FindField("Shape"))
  if (shpfld.IsVisible.Not) then
    shpfld.SetVisible(shpfld.IsVisible.Not)
    WasNotVisible = TRUE
  else
    WasNotVisible =FALSE
  end

  ' This will only export the selected shapes, if any have been selected

  av.ShowMsg("Creating new shape file from"++ t.GetName)
  anFTab = tbl.Export(def, Shape, tbl.GetSelection.Count > 0)
  anFTab.SetEditable(true)

  if (WasNotVisible) then
    shpfld.SetVisible(FALSE)
  end

  ' Find the shape field

  theSField = anFTab.FindField("Shape")

  ' Set up status bar

  av.ShowMsg("Generalizing new shapes from theme"++ t.GetName)
  theSize = anFTab.GetNumRecords

  ' Initialize counters for reporting

  count = 0
  oldCount = 0
  newCount = 0

  ' Process each shape

  for each r in anFTab

    av.SetStatus(100*count/theSize)

    ' Retrieve shape

    if (theClassName = "Polygon") then
      theOShape = Polygon.MakeNull
    else
      theOShape = Polyline.MakeNull
    end
    anFTab.QueryShape(r,theView.GetProjection,theOShape)

    ' Make a list for collecting the new set of point lists

    theNShape = List.Make

    ' Process each list in the ListofListofPoints

    for each p in theOShape.AsList

      ' Make a list for collecting the vertices to keep

      theNList = List.Make

      oldCount = oldCount + p.Count

      ' Set up stack for holding a set of lists of form {index,point}

      theStack = Stack.Make

      ' Add first point to the list being assembled, and make it the 
anchor

      theAnchor = p.Get(0)
      if (hasPrj) then
        theNList.Add(theAnchor.ReturnUnProjected(thePrj))
      else
        theNList.Add(theAnchor)
      end
      aIndex = 0

      ' Add last point to stack

      fIndex = p.Count - 1
      theStack.Push({fIndex,p.Get(fIndex)})

      ' Process the points with the Douglas - Peucker approach

       while (theStack.IsEmpty.Not)

        ' Process from Anchor (beginning) to Floater (top of stack)

        fIndex = theStack.Top.Get(0)

        ' Create graphic shape for comparing distances

        theVect = Line.Make(theAnchor,theStack.Top.Get(1))

        ' Initialize values for comparison

        theMax = theTol
        mxIndex = 0

        ' Find distance

        for each i in (aIndex + 1) .. (fIndex - 1) ' process middle points

          dist = theVect.Distance(p.Get(i))

          ' Check against previous maximum

          if (dist >= theMax) then  ' this point is out of the tolerance
            theMax = dist
            mxIndex = i
          end

        end

        ' If a vertex is found outside the current tolerance corridor,
        ' push it onto the stack

        if (mxIndex > 0) then 'a new floater has been found
          theStack.Push({mxIndex,p.Get(mxIndex)})
        else  ' add floater to the list for the new shape and move forward
          if (hasPrj) then
            theNList.Add(theStack.Top.Get(1).ReturnUnProjected(thePrj))
          else
            theNList.Add(theStack.Top.Get(1))
          end

          ' Make the floater the new anchor

          theAnchor = theStack.Pop.Get(1)
          aIndex = fIndex

        end

      end

      ' Check for collapsed polygons

      if ((theNList.Count < 4) and (theClassName = "Polygon")) then
         ' Revert to original polygon
        if (hasPrj) then
          theNList = List.Make
          for each pt in p
            theNList.Add(pt.ReturnUnProjected(thePrj))
          end
        else
          theNList = p
        end

      end

      ' Finish up individual lists

      newCount = newCount + theNList.Count
      theNShape.Add(theNList)

    end

    ' Finish up the shape

    if (theClassName = "Polygon") then
      theNShape = Polygon.Make(theNShape)
    else
      theNShape = Polyline.Make(theNShape)
    end
    anFTab.SetValue(theSField, r, theNShape)
    count = count + 1

  end

  ' Clear status message

  av.ClearMsg
  av.ClearStatus

  ' Report on counts of old and new theme

  MsgBox.Info("The old theme with"++ oldCount.AsString ++
      "vertices has been generalized to"++
        newCount.AsString++"vertices","Results")

  ' Create a theme and add it to the View

  anFTab.SetEditable(false)
  fthm = FTheme.Make(anFTab)
  theView.AddTheme(fthm)

end

' Bring the View to the front

theView.GetWin.Activate




"Matt Lynch" <matt at terraEngine.com>


"Matt Lynch" <matt at terraEngine.com>
Enviado por: postgis-users-bounces at postgis.refractions.net
28/08/2003 17:45
Por favor, responda a PostGIS Users Discussion

 
        Para:   <postgis-users at postgis.refractions.net>
        cc: 
        Asunto: [postgis-users] How to generalize or simplify a Polygon


Hi,

I have an application where I will be accepting a shape file (polygon)
to define a region of permitted access for a user which will be imported
to postgis.  That shape file has the potential to have many points if
the region defined has a border like a river or lake.  Are there any
functions or tools available that would take a polygon shape and reduce
or generalize the number of points in the polygon?  The resulting
polygon must cover at least as much area as the original.

Thanks,

Matt

---
Outgoing mail is certified Virus Free.
Checked by AVG anti-virus system (http://www.grisoft.com).
Version: 6.0.507 / Virus Database: 304 - Release Date: 8/4/2003
 



_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users