Android Question Polygon around multiple lat/lon points

[RB Smissaert]

I need to write some code to get a polygon (list of lat/lng points) that surrounds a number of lat/lng points on a map.
This polygon needs to surround the points tightly, so a simple rectangle (which would be very easy to code) won't do.
This link looks promising:
https://github.com/mapbox/concaveman
And I wonder if maybe somebody of this group has already coded this for B4A.
In that case I would be very interested to see that code.

RBS

 

[RB Smissaert]

Didn't think it would be a problem but have tested now for the positive to negative longitude change we do have in the UK and that is indeed fine.
Will test now for your mentioned wraparounds.
RBS

 

[emexes]

Will test now for your mentioned wraparounds.

The problem will occur if degrees longitude is eg 358 rather than -2.

I have a vague recollection that most Mod operators don't play well with negative numbers.

Angle = (Angle + 180) Mod 360 - 180

feels like it should work.

Lol I'm just watching an English TV program here, mentioned something about a car worth 70 million.

 

[RB Smissaert]

The problem will occur if degrees longitude is eg 358 rather than -2.

I have a vague recollection that most Mod operators don't play well with negative numbers.

Angle = (Angle + 180) Mod 360 - 180

feels like it should work.

Lol I'm just watching an English TV program here, mentioned something about a car worth 70 million.

> English TV program here, mentioned something about a car worth 70 million
Yeah, just sold that car, shame I didn't realise it was worth that much.

RBS

 

[RB Smissaert]

The problem will occur if degrees longitude is eg 358 rather than -2.

I have a vague recollection that most Mod operators don't play well with negative numbers.

Angle = (Angle + 180) Mod 360 - 180

feels like it should work.

Lol I'm just watching an English TV program here, mentioned something about a car worth 70 million.

I had a go at testing for the situation where longitude changes from eg: 179.99999 to -179.99999 eg when
moving from NE Siberia to Alaska.
I am using OSM maps through this code:

Open Street Map viewer - GPS

Hi, This b4Xlib contains a custom view (cvMap) which can display Open Street Map. The tiles are retrieved from the internet and cached in a database. You can add shapes and images on the map. UI : - Lat/lng Center of the map - Zoom Level - Compass Direction with rotation - Scale - Button...

www.b4x.com

I had to make some alterations to the above code as it gave me the wrong longitude values when moving as
described above.

This is the code in the cvMap class showing may alterations. It will now show the right longitude values, but I am not
sure this is the right approach as more changes are needed as I am unable to draw a polygon covering the above
change in longitude.

'return lng for X tile and offset
Public Sub TileXToLng(aX As Double, aOffsetX As Double, aZoom As Int) As Double

    Dim n As Double = Power(2, aZoom)
    Dim d As Double = ((aX + aOffsetX / MapUtilities.cTileSize) / n) * 360.0
    
    'for some reason working with extreme longitudes, eg 179.99 and -179.99
    'this adjustment is needed to get the right longitude
    'otherwise we only need: d -180
    '----------------------------------------------------------
    If d < 360 Then
        If d < 0 Then
            Dim v As Double = 180 + d
        Else
            Dim v As Double = d -180
        End If
    Else
        Dim v As Double = - 180 + (d - 360)
    End If
    
    'this was the old code as explained above
    '--------------------------------------------------------------------------------
    'Dim v As Double = ((aX + aOffsetX / MapUtilities.cTileSize) / n) * 360.0 - 180.0
    
    'Log("TileXToLng, d: " & d)
    'Log("TileXToLng, MapUtilities.cTileSize: " & MapUtilities.cTileSize)
    'Log("TileXToLng, aOffsetX / MapUtilities.cTileSize: " & (aOffsetX / MapUtilities.cTileSize))
    'Log("TileXToLng, aX: " & aX & ", aOffsetX: " & aOffsetX & ", aZoom: " & aZoom & ", v: " & v)
    
    Return v
    
End Sub

In any case this will need a new thread as the above problem is a very different problem than the problem
of drawing a border around a set of map points.
Will post this new thread once I have a bit better idea what is going on.

RBS

 

[emexes]

I think the trick would be to trust that the point cloud will be at-most 180 degrees wide, and thus only include (at-most) one of the two obvious wraparounds ie either meridian 0 or meridian 180.

If point cloud includes meridian 0, then use longitudes -180 to 179.999
Angle = (Angle + 180) Mod 360 - 180

If point cloud includes meridian 180, then use longitudes 0 to 359.999
Angle = (Angle + 360) Mod 360

If point cloud doesn't include either meridian, then doesn't matter which longitude scheme is used

What could possibly go rwong?

 

[emexes]

I figure the widest point cloud would be Russia, from 19° to 191°, which crosses the meridian at 180°

I realised that point clouds > 180° still work ok, as long as they only cross one of the 0° and 180° meridians.

So it'd be a case of sorting all the longitudes (in 0 to 359.999° form) and then finding the largest longitudinal gap, which will be either the largest gap between two consecutive longitudes, or the gap between the last longitude and the first longitude. If the largest gap is the latter, ie across the 0° meridian, then keep using the 0 to 359.999° form, otherwise use the -180 to 179.999° form.

Lol I like that you're already eyeing worldwide domination in whatever it is that your app does.

 

[emexes]

If the largest gap is the latter

Ugh that works great as long as the point cloud doesn't cross both meridians.

I think the solution is to temporarily offset the longitudes by the middle (or any) longitude within the largest longitudinal gap, which will shift the point cloud to not cross the wraparound, then do the bounding polygon, then remove the offset from the longitudes.

 

[RB Smissaert]

I think the trick would be to trust that the point cloud will be at-most 180 degrees wide, and thus only include (at-most) one of the two obvious wraparounds ie either meridian 0 or meridian 180.

If point cloud includes meridian 0, then use longitudes -180 to 179.999
Angle = (Angle + 180) Mod 360 - 180

If point cloud includes meridian 180, then use longitudes 0 to 359.999
Angle = (Angle + 360) Mod 360

If point cloud doesn't include either meridian, then doesn't matter which longitude scheme is used

What could possibly go rwong?

They are all good points and will bear those in mind when I come to testing ConvexHull at those extreme longitude values.
Problem is that I am not able to test easily yet as I have some troubles with my map code at those extreme values.
I test by drawing a polygon on the map, fill that polygon with random map points and then run ConvexHull on those points.
Between the longitude values of about -179.999 to -175 (even when not crossing the +/- line) there is a problem drawing that
polygon as the map code as doesn't pick up the right drawing points (selected by touching the map or clicking a button when
the map center is at the required point).

I think the problem might be in this bit of map code:

'-----------------------------------------------------------------------------------------------------------------
'this can give negative values for TMapPoint.fX if providing eg -179.999 for longitude and that doesn't seem right
'-----------------------------------------------------------------------------------------------------------------
'return  TMapPoint from TMapLatLng
Public Sub LatLngToPoint(aLatLng As TMapLatLng) As TMapPoint
    
    'https://www.netzwolf.info/osm/tilebrowser.html?marker.x=12&marker.y=189&tx=1&ty=0&tz=1&ts=256#tile
    Dim v As B4XView = PointToTile(0, 0)
    Dim txy As TMapTileXY = v.tag
    Dim gpx As Double = (txy.fX + (-v.Left / MapUtilities.cTileSize)) / Power(2, fMap.fZoomLevel)
    Dim gpy As Double = (txy.fy + (-v.top / MapUtilities.cTileSize)) / Power(2, fMap.fZoomLevel)
    Dim tx As TMapTileNumberOffset = LngToTileXPlusOffset(aLatLng.fLng, fMap.fZoomLevel)
    Log("cvMap, LatLngToPoint, tx.fTile : " & tx.fTile & ", tx.fOffset : " & tx.fOffset)
    Dim ty As TMapTileNumberOffset = LatToTileYPlusOffset(aLatLng.fLat, fMap.fZoomLevel)
    Dim ppx As Double = (tx.fTile + (tx.fOffset/MapUtilities.cTileSize)) / Power(2, fMap.fZoomLevel)
    Log("cvMap, LatLngToPoint, ppx : " & ppx)
    Dim ppy As Double = (ty.fTile + (ty.fOffset/MapUtilities.cTileSize)) / Power(2, fMap.fZoomLevel)
    Dim p As TMapPoint = MapUtilities.initPoint(fTilesCount * MapUtilities.cTileSize * (ppx - gpx), fTilesCount*MapUtilities.cTileSize * (ppy - gpy))
    
    Return p
    
End Sub

But as said, this is really a different problem than the topic of this thread, so I will have to start a new thread for that if I can't fix this problem.

RBS

 

[MasterGy]

I have a suggestion.

Given the raw, original longitude/latitude points.
Find the limits.
Xmin, Xmax, Ymin, Ymax.

Find the middle (average) of the limits
refX = (Xmin+Xmax)/2
refY = (Ymin+Ymax)/2

And load the points into the polygon generator like this:

newX = difang(refX,originalX)
newY = difang(refY,originalY)

Sub  difang (angle1 As Float, angle2 As Float) As Float
    Dim pi as Float = 3.1415926 'if you use radian 0-pi
    Dim pi as Float = 180 'if you use degree
  

    Dim diff As Float
    Dim c As Float =  2 * pi
    Do Until angle1 < c:        angle1 = angle1 - c:    Loop
    Do Until angle1 >= 0:        angle1 = angle1 + c:    Loop
    Do Until angle2 < c:        angle2 = angle2 - c:    Loop
    Do Until angle2 >= 0:        angle2 = angle2 + c:    Loop
    diff = angle1 - angle2
    If diff < -pi Then diff = diff + 2 * pi
    If diff > pi Then diff = diff - 2 * pi
    Return  diff
End Sub

You create a reference point within your points. That will be 0. The 'difang' sub calculates the difference between two degrees, and takes into account that 0=360. For example, difang(350,10) will not be 340, but -20. This ensures that the difference between -179 and 179 degrees is not 358 degrees, but 2. I can't test it, but I would try in that direction.

This might solve the problem.

 

[RB Smissaert]

I have a suggestion.

Given the raw, original longitude/latitude points.
Find the limits.
Xmin, Xmax, Ymin, Ymax.

Find the middle (average) of the limits
refX = (Xmin+Xmax)/2
refY = (Ymin+Ymax)/2

And load the points into the polygon generator like this:

newX = difang(refX,originalX)
newY = difang(refY,originalY)

Sub  difang (angle1 As Float, angle2 As Float) As Float
    Dim pi as Float = 3.1415926 'if you use radian 0-pi
    Dim pi as Float = 180 'if you use degree
 

    Dim diff As Float
    Dim c As Float =  2 * pi
    Do Until angle1 < c:        angle1 = angle1 - c:    Loop
    Do Until angle1 >= 0:        angle1 = angle1 + c:    Loop
    Do Until angle2 < c:        angle2 = angle2 - c:    Loop
    Do Until angle2 >= 0:        angle2 = angle2 + c:    Loop
    diff = angle1 - angle2
    If diff < -pi Then diff = diff + 2 * pi
    If diff > pi Then diff = diff - 2 * pi
    Return  diff
End Sub

You create a reference point within your points. That will be 0. The 'difang' sub calculates the difference between two degrees, and takes into account that 0=360. For example, difang(350,10) will not be 340, but -20. This ensures that the difference between -179 and 179 degrees is not 358 degrees, but 2. I can't test it, but I would try in that direction.

This might solve the problem.

Thanks, will look at that when I am able to test on the map.
RBS

 

[emexes]

Given the raw, original longitude/latitude points.
Find the limits.
Xmin, Xmax, Ymin, Ymax.

This is certainly a clearer and easier-to-understand explanation than my hand-wavy attempt, but:

1/ no need to do Y axis (latitude) because it doesn't wrap around like longitude does eg when you cross over +90° latitude at north pole, you are not magically transported to -90° at south pole.

2/ finding the longitude limits still fails if the point cloud includes the wraparound median eg London ±1° east/west (approx. ±55 miles) would, if using 0..359.999 notation, be 359° to 1° and the middle of that range be 180° ie the other side of the world. Using the "largest gap" method of post #87 would work correctly.