Point, Line, PolyLine, Polygon

[Paul Higgins]

Yes, it works great. Thank you.
Question on the : intersectionPoint( LS2, 2 ) what does the ",2" do? If I delete it it gives same result.

Also, is there some documentation to explain what all the methods do?

[Jan-Erik Lärka]

The method usually loop over the points defined for the shape from first to last.
The added 2 tell it to not check the first line segment between the first and second point, but rather start at the second line segment.
The second argument is useful if you want to find the possible second, third, ... intersection point as the code just return the first one found as it is.

I haven't begun to write documentation for it, as I daily add and adjust more functionality to code that this is based on, at work, as I have to handle more and more all the time. This is my attempt to clean up the code and structure it better.

"intersects" return .true or .false if they intersect
"intersect" return the line segments that intersect as a number dot number
"intersectionPoint" now return the point of intersection

Note! The PolyShape class is unaware of the curvature of the earth, so it is inaccurate for larger distances. PolyShape calculate the intersection point for a flat surface.

Polygon.cls can contain and handle for example holes in the shape

WGS84.cls handle functionality for earth curvature suitable with google kml files and gps. Need help to find code here to calculate intersection point for a curved surface such as the earth.

---

Feel free to test each class and method, document what they can and can't do.
I've tested a package that "document" classes automatically, and it seem to do a nice job, but one need to add useful insights, examples and explain what each class should be used for.

I would also be pleased if for example method "offset" could be completed to return shapes with
a) correct offset for ends when used with both side offset, not only inside or outside.
b) clean up of line segments that overlap due to large offset distance.

[Paul Higgins]

Jan, Thank you. I have been using OORexx for nearly 40 years and I am getting better all the time but when I look at your code I am still an amateur. Keep up the great work.

Cheers,

[Jan-Erik Lärka]

Thanks,

I'm new to ooRexx as it wasn't ported to OS/2 and Object Rexx wasn't the default interpreter on OS/2 so one could not rely on that is was enabled, while classic rexx was.

May not do things as efficient as one could, but that is because how and what the descriptions/documentation say how to write code.

Polygon and related code/classes mostly use techniques found in school book examples though there usually written in other languages.

I think that "rexx" would be useful in education as at least pseudo code, but often usable code to create understanding.

[Paul Higgins]

Jan, if you have time look at this.

I define two polygons: see the attachment, the first polygon is in blue the second in yellow. They have two intersection points.

curr_poly =.PolyShape~new('117.519585,8.643308,0       117.522866,8.583716,0   117.590189,8.587334,0   117.586919,8.646926,0   117.519585,8.643308,0')
big_poly = .PolyShape~new('117.526397,8.600475,0        117.47546,8.500079,0    117.535704,8.470178,0   117.586657,8.570574,0   117.526397,8.600475,0')
say 'intersect='curr_poly~intersect(big_poly)
IP0 = curr_poly~intersectionPoint(big_poly)
say 'IP0 intersectionPoint='IP0~y IP0~x 
IP0 = curr_poly~intersectionPoint(big_poly,2)
say 'IP0 intersectionPoint='IP0~y IP0~x 

The output is:
intersect=1.1
IP0 intersectionPoint=8.59255535 117.522379
IP0 intersectionPoint=8.59255535 117.522379

How do I see the 2nd intersection point ?
Thanks

[Jan-Erik Lärka]

Hadn't removed a thing in delta...
well, here's some more to play with
2RectPolyto test PolyShape.cls
odc WGS84.cls to list classes

[Paul Higgins]

Jan,

Yes that fixed it thanks.

BTW I changed the self~decimals = 3 to 6 in PolyShape.Cls so the tostring gives the full 6 decimal places. Having 3 decimal places is not too useful for navigation . I think you should make it 6 as most users would not know how to change it.

Will checkout your other files later.

[Jan-Erik Lärka]

It's deliberate, as I want you to use:
a) Polygon (and WGS84 that already use more decimals).
b) LS0~decimals = 6

[Paul Higgins]

Jan-Erik,
I seem to still have a problem finding the intersectionpoint of two lines when the angle between them is small. Here is what I found:

BG_LN =.PolyShape~new('-77.619237,20.631839 -77.619237,20.698505') 
say 'BG_LN='BG_LN
LG_LN =.PolyShape~new('-77.617041,20.698442 -77.621432, 20.631902')  
say 'LG_LN='LG_LN  
say 'intersect='LG_LN~intersect(BG_LN)                  -- 
say 'intersects='LG_LN~intersects(BG_LN)  
IP1 = LG_LN~intersectionPoint(BG_LN)
if .Nil <> IP1 then  
   say 'intersectionPoint='IP1~y IP1~x
else 
   say 'intersectionPoint is .Nil'

the result is:
BG_LN=a Line
LG_LN=a Line
intersect=1.1
intersects=1
intersectionPoint is .Nil

if I decrease the angle of the BG_LN it works
LN =.PolyShape~new('-77.619237,20.631839 -77.622473,20.698414')
BG_LN=a Line
LG_LN=a Line
intersect=1.1
intersects=1
intersectionPoint=20.6510311 -77.6201699

[Paul Higgins]

OK I found the problem: at the end of the intersectionpoint method you do compares and some fail because they are are not formatted to the same number of decimal places. Fixed it by:

 Do j = 1 To shape~point~Items - 1
 (shape~point[ j ]~x)~format(4,6)~strip 
 (shape~point[ j +1 ]~x)~format(4,6)~strip 
 (shape~point[ j ]~y)~format(3,6)~strip 
 (shape~point[ j +1 ]~x)~format(3,6)~strip 
 ...
 px = ( d2~x * d1~z - d1~x * d2~z ) / d12
 py = ( d1~y * d2~z - d2~y * d1~z ) / d12
 px = px~format(4,6)~strip
 py = py~format(3,6)~strip

 Attached new polyshape.cls with fix

[Paul Higgins]

I changed all formats in my previous post to:
~format(18,self~decimals)~strip
to make it compatible with all decimal values.

But to make that test above work I had to change the self~decimals = 6 for the Polyshape class and use .Polyshape for both lines because if I use WGS84 in my test for both the lines it always returns .Nil
like this:

BG_LN =.WGS84~new('-77.619237,20.631839 -77.619237,20.698505')
say 'BG_LN='BG_LN
LG_LN =.PolyShape~new('-77.617041,20.698442 -77.621432, 20.631902')
say 'LG_LN='LG_LN  
say 'intersect='LG_LN~intersect(BG_LN)                  -- 
say 'intersects='LG_LN~intersects(BG_LN)  
IP1 = LG_LN~intersectionPoint(BG_LN)
if .Nil <> IP1 then  
   say 'intersectionPoint='IP1~y IP1~x
else 
   say 'intersectionPoint is .Nil'

returns:
BG_LN=a Line
LG_LN=a Line
intersect=0
intersects=0
intersectionPoint is .Nil

So I can't use WGS84 or I am I missing the point?

[Jan-Erik Lärka]

/**/

BG_LN =.PolyShape~new
BG_LN~decimals = 12 --set attribute that affect preparation/interpretation.
BG_LN~prepare('-77.619237,20.631839 -77.619237,20.698505') 
say 'BG_LN='BG_LN
LG_LN =.PolyShape~new
LG_LN~decimals = 14 -- set attribute to different value for this.
LG_LN~prepare('-77.617041,20.698442 -77.621432, 20.631902')  
say 'LG_LN='LG_LN  
say 'intersect='LG_LN~intersect(BG_LN)                  -- 
say 'intersects='LG_LN~intersects(BG_LN)  
IP1 = LG_LN~intersectionPoint(BG_LN,, .True)
if .Nil <> IP1 then  
   say 'intersectionPoint='IP1~y IP1~x
else 
   say 'intersectionPoint is .Nil'

::requires 'PolyShape.cls'

[Jan-Erik Lärka]

Added use of self~format( ... ) to calculation for px and py instead. I should have thought about it as the method crop that I develop for work also require it.

[Paul Higgins]

Thanks for the update. But I still can not use WGS84.

[Jan-Erik Lärka]

/* PolyShape can't handle WGS84 shapes as WGS84 is based on PolyShape and not the other way around, with 'an Array' in between */

BG_LN =.wgs84~new('-77.619237,20.631839 -77.619237,20.698505') 
say 'BG_LN='BG_LN~toString( ,,,'' ) --OBS!! Compare with LG_LN
LG_LN =.wgs84~new
LG_LN~decimals = 6 --Set this as well if you replace wgs84 with polyshape on the line above
LG_LN~prepare('-77.617041,20.698442 -77.621432, 20.631902')  
say 'LG_LN='LG_LN~toString( ,,.True )  --OBS!! Use toString to get formatted output
say 'intersect='LG_LN~intersect(BG_LN) --Set to BG_LN~shapes[1] if BG_LN = .wGS84 & LG_LN = .PolyShape
say 'intersects='LG_LN~intersects(BG_LN)  --Set to BG_LN~shapes[1] if BG_LN = .wGS84 & LG_LN = .PolyShape
IP1 = LG_LN~intersectionPoint(BG_LN) --Set to BG_LN~shapes[1] if BG_LN = .wGS84 & LG_LN = .PolyShape
if .Nil <> IP1 then
say 'intersectionPoint='ip1~toString( .false ) --See toString in .Point, .PolyShape etc.
else 
   say 'intersectionPoint is .Nil'

::requires 'WGS84.cls'
::requires 'PolyShape.cls'

[Paul Higgins]

Jan-Erik,
Something is really confusing to me. Why do you reverse the x and y for a point. For example:

pt = .point~new("117.501452,8.490863")
say 'PT='pt~x pt~y

the result is:
PT=8.490863 117.501452

[Jan-Erik Lärka]

Google use "Lat,Lng" for an individual point, that is in reversed order, while not for polylines and polygons.
longitude = x
latitude = y
elevation = z or e
<...> indicate optional parameter

The point class can use/handle/interpret these formats:
pnt = .point~new( x, y<, z> ) -- each set individually
pnt2 = .point~new( '(x,%20y<,%20z>)' ) -- interpret string as point found in web pages
pnt3 = .point~new( '(x,y<,z>)' ) -- interpret string as point
pnt4 = .point~new( 'y,x<,e>' ) -- interpret point from google string seen in google maps etc.

So use paranthesis around, such as '(x,y)', to interpret in longitude,latitude order.