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
Find the union of two (2) polygons
-
- Posts: 17
- Joined: Jan 08, 2016 20:33
Re: Find the union of two (2) polygons
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!
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!
Re: Find the union of two (2) polygons
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
-
- Posts: 17
- Joined: Jan 08, 2016 20:33
Re: Find the union of two (2) polygons
Thanks. I’ll have a look at this when I get a chance.
Your help is greatly appreciated.
Regards
Andrew
Your help is greatly appreciated.
Regards
Andrew