View previous topic :: View next topic |
Author |
Message |
snipsnip
Joined: 10 Feb 2014 Posts: 74 Location: Melbourne
|
Posted: Sat Nov 30, 2019 12:00 pm Post subject: point in polygon function and kml placemark test |
|
|
This is a work in progress, it works but the poly point array must start and end with the same coords.
I'm having a hard time getting this to happen in the next loop. ie closing the polygon
taken from one of the many pnp ray tracing algorithms on line
as it is now, its enough for what i need.
The test is just a static one - you'll need to capture the serial output and send it to a kml file
YAT (yet another terminal) works well for this
Code: | Function PIP(byval X as single , byval Y as single) as byte
' index is the polygon point array count
' fencelon() and fencelatt() global arrays of the x,y coords - fencelon is the "Y" axis
' x,y are the test point coords - x=latt y=lon
local f1 as single, f2 as single,f3 as single,f4 as single,f5 as single ,f6 as single
local cn as byte
local i as byte, j as byte ,q as byte,r as byte
i=1 : cn = 0
index = index-1
j=index
for i = 1 to index
q=0 : r=0
if fencelon(i)< Y And fencelon(J)>= Y then q = 1
if fencelon(J)< Y And fencelon(i)>= Y then r=1
If q=1 or r=1 then
f1 = Y - fencelon(i) : f2 = fencelon(j) - fencelon(i) : f3 = fencelatt(j) - fencelatt(i)
f4 = f1/f2
f5 = f4 * f3
f6 = fencelatt(i) + f5
if f6 < X Then cn = cn+1
End If
j=i
next i
index = index+1
cn = cn mod 2
if cn > 0 then cn =1 ' even numbers are outside the polygon
PIP = cn
End Function
'-----------------------------------------------------------------------------------------------------------------------------------------
'-----------------------------------------------------------------------------------------------------------------------------------------
sub PIPtestKML(byval OffSet as single, byval resolution as single)
'decimal place / grid resolution / offset
'2 0.01 1.1132 km
'3 0.001 111.32 m
'4 0.0001 11.132 m
'5 0.00001 1.1132 m
'go easy with these - easy to create many millions of points, google earth will have a hissy fit.
' check the terminal program for max number of displayed lines.
'have been using YAT terminal for this test, - need more than 80 columns and can change the buffered displayed lines easy.
'find max lat and min latt, max lon min lon from stored fencelatt,lon array and add (offset) to them
'this test wont work accross the equator or the prime meridian
local n as byte ,m as byte
local maxlon as single, minlon as single, maxlatt as single, minlatt as single , original as single
'set the focus point for initial Earth view and kml header stuff
print "<?xml version=";chr(34);"1.0";chr(34);" encoding=";chr(34);"UTF-8";chr(34);"?>"
print "<kml xmlns=";chr(34);"http://www.opengis.net/kml/2.2";chr(34);
print " xmlns:gx=";chr(34);"http://www.google.com/kml/ext/2.2";chr(34);">"
print "<Document>"
print "<LookAt> <longitude>";
print fencelon(1);
print "</longitude> <latitude>";
print fencelatt(1);
print "</latitude> <range>1300.000000</range> </LookAt>"
print "<Folder>"
''' gegerate the FENCEPOINTS / POLYPOINTS
for n = 1 to index
print "<Placemark>"
print "<name>";"FP";n;"</name>"
print"<Point> <coordinates>"
print fencelon(n);
print ",";
print fencelatt(n);
print "</coordinates></Point></Placemark> "
next n
'find max and mins then add the offset value to generate the test box
maxlatt = fencelatt(1) : minlatt = fencelatt(1): maxlon = fencelon(1) : minlon = fencelon(1)
for n = 1 to index
if maxlatt <= fencelatt(n) then maxlatt = fencelatt(n)
if minlatt >= fencelatt(n) then minlatt = fencelatt(n)
if maxlon <= fencelon(n) then maxlon = fencelon(n)
if minlon >= fencelon(n) then minlon = fencelon(n)
next n
minlatt = minlatt - offset
maxlatt = maxlatt + offset
minlon = minlon - offset
maxlon = maxlon + offset
''' run the test grid through the point in polygon(pip) function and generate placemark pins
'' original = minlatt
while minlatt < maxlatt and minlon < maxlon 'this sets the test boundry box - it's just easier to do a test square
original = minlatt
while minlatt < maxlatt ' print the x axis
m = pip(minlatt, minlon)
if m< 1 then ' just display the outside points
print "<Placemark>"
print "<name>";m;"</name>"
print" <Point><coordinates>"
print minlon;
print ",";
print minlatt;
print "</coordinates></Point></Placemark>"
end if
minlatt = minlatt + resolution
wend
minlatt = original
minlon = minlon + resolution
wend
print "</Folder></Document></kml>"
do
loop ' stop here - need time to save to file in terminal
end sub
sub preload()
fencelatt(1)=-37.702547
fencelatt(2)=-37.702413
fencelatt(3)=-37.701884
fencelatt(4)=-37.702027
fencelatt(5)=-37.702547 'start repeat
fencelon(1)=147.840474
fencelon(2)=147.839278
fencelon(3)=147.839317
fencelon(4)=147.840353
fencelon(5)=147.840474 'start repeat
index = 5
end sub
|
call PIPtestKML(0.0005,0.00005) generated the kml below
with the preload sub lat lon values - for an idea of resolution
|
|
Back to top |
|
|
albertsm
Joined: 09 Apr 2004 Posts: 5913 Location: Holland
|
Posted: Thu Dec 05, 2019 9:59 am Post subject: |
|
|
thank you for sharing. it looks great, i needed some time to understand what it is
i use kml files too when i walk to make a log file which i later look at with google earth.
I also did some tracking apps with gsm/wifi.
But it took me some time to understand the purpose. this is not just a square but a polygon.
great for fencing.
thanks! _________________ Mark |
|
Back to top |
|
|
Matrixx
Joined: 30 Nov 2005 Posts: 158
|
Posted: Wed Aug 19, 2020 6:21 pm Post subject: |
|
|
Than you for publish this great example! : )
I was in the need of doing the same polygon test for one of my projects.
However, reading in detail the author page, I noticed the I & J must evaluate all the vertex.
I compared the original Author code in C++, also the original code of the same Author in Fortran.
in the function Code: | Function PIP(byval X as single , byval Y as single) as byte |
You have
Code: | index = index-1
J = index |
when the function starts, so this leave last vertex out of the loop.
Maybe this is the reason you have issues about need to close the polygon (repeat the first vertex);
in fact the function do not need to repeat the first vertex, once it needs to cycle thru all the vertex.
The vertex order CW or CCW does not matter.
i.e. having a rectangle, you need to pass only 4 vertex to the function (index = 4) and not repeat the first one (index = 5).
So my 0.5 cents are a couple of things need to change:
1. do the loop from 1 to index (do not subtract 1 from index )
2. assign the value of J = N at the beginning, before the loop starts.
This way, the first loop iteration will evaluate I=1, J=N so effectively evaluating first and last vertex.
The rest of the loop iterations evaluates what else is in between.
It is working great.
Example: testing rectangle with 4 vertex:
|
|
Back to top |
|
|
|
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum You cannot attach files in this forum You cannot download files in this forum
|
|