Find the union of two (2) polygons

General FreeBASIC programming questions.
Post Reply
Andrew Lindsay
Posts: 17
Joined: Jan 08, 2016 20:33

Find the union of two (2) polygons

Post by Andrew Lindsay »

OK, so after much very much appreciated help with my last conundrum, I now stumble onto the next phase of my journey.

I would like to know how to create a combined outline of two overlapping datasets. I do not want to find the smallest bounding box, but I actually want to merge or join a union of the two datasets to find the resultant single shape.

The datasets will be in simple XY sets, forming a complete closed shape. Each shape, one dataset. I would like to then combine the tow overlapping outlines to form a single outline.

Any assistance with this next challenge would be greatly appreciated.

Regards
Andrew
caseih
Posts: 2157
Joined: Feb 26, 2007 5:32

Re: Find the union of two (2) polygons

Post by caseih »

I did a quick google search and found various algorithms and libraries that could do boolean operations on polygons which you might implement. Nothing existing in FB of course. There are various libraries including GEOS, which has a plain C api that could be called from FB. Here's a union-only algorithm described (but not coded): https://stackoverflow.com/questions/266 ... x-polygons.

This operation is pretty trivial to do in a scripting language I typically use, using the shapely library.

However, maybe coding in FB or any programming language is the wrong approach. Sometimes it's better to use a dedicated tool. Sounds like you are doing GIS analysis, which is what QGIS is meant for. And you can visualize the polygon on a map as well. QGIS is the tool I would use for this rather than code up an algorithm or even writing a script and using a library. I would import your CSVs into QGIS. I would specify either the correct UTM grid, or use no coordinate system at all (assuming your csvs are all consistent in using the same coordinate system). It's trivial in QGIS to merge two or more polygons into a new one. Then I would export the polygon to a csv file. If I needed to do this many times it is possible to make macros in QGIS. QGIS is very powerful but has a steep learning curve. However if you're doing GIS analysis regularly, it's worth your time to learn. I know this likely isn't the answer you are looking for!
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Find the union of two (2) polygons

Post by dodicat »

For fun

Code: Select all


'#cmdline "-exx"  'OK
Const u32max = &hFFFFFFFFul             ' ULong maximum value
#define rnd64u (culngint(rnd*u32max)+(culngint(rnd*u32max) shl 32))
#define rng32(f,l)  clng((culngint(rnd64u) mod (((l)-(f))+1)) + (f))


Type Point
      As Double x,y,z,k
      Declare Operator Cast() As String
End Type

Operator point.cast() As String 'unused
Return "("+Str(x)+","+Str(y)+")"
End Operator

Operator =(a As Point,b As Point) As boolean
Return a.x=b.x And a.y=b.y
End Operator

Function inpolygon(p1() As Point,Byval p2 As Point) As Integer
      #define Winder(L1,L2,p) ((L1.x-L2.x)*(p.y-L2.y)-(p.x-L2.x)*(L1.y-L2.y))
      Dim As Long index,nextindex,k=Ubound(p1)+1,wn
      For n As Long=1 To Ubound(p1)
            index=n Mod k:nextindex=(n+1) Mod k
            If nextindex=0 Then nextindex=1
            If p1(index).y<=p2.y Then
                  If p1(nextindex).y>p2.y Andalso  Winder(p1(index),p1(nextindex),p2)>0 Then wn+=1
            Else
                  If p1(nextindex).y<=p2.y Andalso Winder(p1(index),p1(nextindex),p2)<0 Then wn-=1
            End If
      Next n
      Return wn
End Function

Sub drawpolygon(p() As Point,bdr As Ulong)
      Dim k As Long=Ubound(p)+1
      Dim As Long index,nextindex
      For n As Long=Lbound(p) To Ubound(p)
            index=n Mod k:nextindex=(n+1) Mod k
            If nextindex=0 Then nextindex=Lbound(p)
            Line (p(index).x,p(index).y)-(p(nextindex).x,p(nextindex).y),bdr
      Next
End Sub


#macro cleanup(a)
Scope
      Redim As Typeof(a) b(Lbound(a) To Ubound(a))
      b(Lbound(b))=a(Lbound(a))
      Dim As Long flag,count=1
      For n1 As Long=Lbound(a)+1 To Ubound(a)
            flag=0
            For n2 As Long=Lbound(a) To  n1-1
                  If a(n1)=a(n2) Then flag=1:Exit For
            Next n2
            If flag=0 Then
                  b(Lbound(b)+count)=a(n1)
                  count+=1
            End If
      Next n1
      Redim Preserve b(Lbound(a) To count+Lbound(a)-1)
      Redim a(Lbound(b) To Ubound(b))
      For n As Long=Lbound(b) To Ubound(b)
            a(n)=b(n)
      Next
End Scope
#endmacro


#macro IsInArray( N,s,b)
Do
      For i As Long = Lbound(s) To Ubound(s)
            If n=s(i) Then b=true:Exit Do
      Next i
      b=false
      Exit Do
Loop
#endmacro

#macro AddToArray(n,s)
Do
      Dim As boolean b
      IsInArray( n,s,b)
      If b =true Then  Exit Do
      Redim Preserve s(Lbound(s) To Ubound(s)+1)
      s(Ubound(s))=n
      Exit Do
Loop
#endmacro

#macro ArraysUnion(a,b,u)
Scope
      Dim As Integer ub=Ubound(a)
      Redim u(Lbound(a) To ub)
      For n As Long=Lbound(a) To ub
            u(n)=a(n)
      Next
      Redim Preserve u(Lbound(a) To ub+(Ubound(b)-Lbound(b)+1))
      For n As Long=Lbound(b) To Ubound(b)
            u(ub+n+1)=b(n)
      Next n
      cleanup(u)
End Scope
#endmacro


#macro IntersectArrays(s,T,U)
Do
      Dim As boolean b
      For i As Integer = Lbound(s) To Ubound(s)
            Var k=s(i)
            IsInArray(k,t,b)
            If b=true Then : AddToArray(k,U):End If
      Next i
      Exit Do
Loop
#endmacro
'======================
Dim As Point poly1(1 To 5)={(20,20),(400,50),(480,400),(200,550),(30,300)}
Dim As Point poly2(1 To 4)={(410,20),(700,50),(580,400),(300,550)}',(30,300)}
'array 1 setup
Dim As Long lim=300*25,count=-1
Redim As Point a1(lim)

For m As Long=1 To lim
      Var x=rng32(0,500),y=Int(Rnd*600)
      If inpolygon(poly1(),Type(x,y)) Then
            count+=1
            a1(count)=Type<Point>(x,y)
      End If
Next m
Redim Preserve a1(count)
cleanup(a1)

'array 2 set up
count=-1
lim=500*25
Redim As Point a2(lim)
For m As Long=1 To lim
      Var x=rng32(300,800),y=Int(Rnd*600)
      If inpolygon(poly2(),Type(x,y)) Then
            count+=1
            a2(count)=Type<Point>(x,y)
      End If
Next m
Redim Preserve a2(count)
cleanup(a2)
'=========================
Print "please wait"
Redim As Point un()
Redim As Point isect()

Screen 19,32,,64
drawpolygon(poly1(),Rgba(255,255,255,200))
drawpolygon(poly2(),Rgba(255,255,255,200))
Color Rgba(255,0,255,255)
Print "array 1  "+Str(Ubound(a1)+1)
For n As Long=Lbound(a1) To Ubound(a1)
      Pset(a1(n).x,a1(n).y),Rgba(255,0,255,255)
Next

Color Rgba(0,255,255,255)
Print "array 2  "+Str(Ubound(a2)+1)
For n As Long=Lbound(a2) To Ubound(a2)
      Pset(a2(n).x,a2(n).y),Rgba(0,255,255,255)
Next

ArraysUnion(a1,a2,un)
Color Rgba(200,200,200,80)
Print "union    "+Str(Ubound(un)+1)
For n As Long=Lbound(un) To Ubound(un)
      Circle(un(n).x,un(n).y),3,Rgba(200,200,200,80)
Next

IntersectArrays(a1, a2,isect)
Color Rgba(255,255,0,100)
Print "intersection "+Str(Ubound(isect)+1)
For n As Long=Lbound(isect) To Ubound(isect)
      Circle(isect(n).x,isect(n).y),5,Rgb(255,255,0)',,,,f
Next

For n As Long=Lbound(isect) To Ubound(isect)
      Draw String(20,60+n*16),Str(isect(n))
      ' print isect(n)
Next

Color Rgb(255,255,255)

Draw String(20,560), "done"
Sleep

 
Andrew Lindsay
Posts: 17
Joined: Jan 08, 2016 20:33

Re: Find the union of two (2) polygons

Post by Andrew Lindsay »

Thanks. I’ll have a look at this when I get a chance.

Your help is greatly appreciated.
Regards
Andrew
Post Reply