golang / geo

S2 geometry library in Go
Apache License 2.0
1.68k stars 182 forks source link

Loop.Intersects gives incorrect results #51

Closed outsidebedisdoor closed 5 years ago

outsidebedisdoor commented 5 years ago

I have two geojson objects Polygon 1: [[[-74.93335056304832,39.283340980346445], [-74.93335056305031,40.30000217770297], [-75.91667795181374,40.30000217770097], [-75.91667795181175,40.28334098034444], [-74.93335056304832,39.283340980346445]]] Polygon 2: [[[-74.3122237229173, 39.9541747627284], [-74.1955570274039, 40.0833414585103], [-74.5497238128233, 40.1958414684098], [-74.7163905245327, 40.0583414352459], [-74.8663905578957, 39.8208413826612], [-74.4163904122584,39.7708413884741], [-74.3122237229173, 39.9541747627284]]]

I converted [][][]float64 array to Loop using Function: func convertFloatListToPolygon(floatList [][][]float64) *s2.Loop { points := floatList[0] pointList := make([]s2.Point, 0) for _, point := range points { s2Point := s2.PointFromLatLng( s2.LatLng{ Lat:s1.Angle(point[1]), Lng:s1.Angle(point[0]), }) pointList = append(pointList, s2Point) } return s2.LoopFromPoints(pointList) }`

poly = convertFloatListToPolygon(p1) comparedPoly = convertFloatListToPolygon(p2) poly.Intersects(comparedPoly)

The results is true. However, when I draw those two polygon in geojson.io, it's not intersected.

Could someone tell me if I use function wrong or there is a bug in intersection?

rsned commented 5 years ago

Both of the loops are invalid because they have duplicate vertices. The duplicate vertex leads to undefined results in most of the resulting calculations.

l1.Validate() = edge 4 is degenerate (duplicate vertex) l2.Validate() = edge 6 is degenerate (duplicate vertex)

On Tue, May 21, 2019 at 5:12 PM Yanhe Chen notifications@github.com wrote:

I have two geojson objects Polygon 1: [[[-74.93335056304832,39.283340980346445], [-74.93335056305031,40.30000217770297], [-75.91667795181374,40.30000217770097], [-75.91667795181175,40.28334098034444], [-74.93335056304832,39.283340980346445]]] Polygon 2: [[[-74.3122237229173, 39.9541747627284], [-74.1955570274039, 40.0833414585103], [-74.5497238128233, 40.1958414684098], [-74.7163905245327, 40.0583414352459], [-74.8663905578957, 39.8208413826612], [-74.4163904122584,39.7708413884741], [-74.3122237229173, 39.9541747627284]]]

I converted [][][]float64 array to Loop using Function:

func convertFloatListToPolygon(floatList [][][]float64) *s2.Loop { points := floatList[0] pointList := make([]s2.Point, 0) for _, point := range points { s2Point := s2.PointFromLatLng( s2.LatLng{ Lat:s1.Angle(point[1]), Lng:s1.Angle(point[0]), }) pointList = append(pointList, s2Point) } return s2.LoopFromPoints(pointList) }

poly = convertFloatListToPolygon(p1) comparedPoly = convertFloatListToPolygon(p2) poly.Intersects(comparedPoly)

The results is true. However, when I draw those two polygon in geojson.io, it's not intersected.

Could someone tell me if I use function wrong or there is a bug in intersection?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/golang/geo/issues/51?email_source=notifications&email_token=ACBHGBEZHP2M6TZQIODQCZ3PWSFVPA5CNFSM4HOPUNO2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4GVCZVFQ, or mute the thread https://github.com/notifications/unsubscribe-auth/ACBHGBEQCEWETWBF3WBHRALPWSFVPANCNFSM4HOPUNOQ .

outsidebedisdoor commented 5 years ago

@rsned Thanks for answering. Is there an efficient way to make a valid loop then?

rsned commented 5 years ago

Looking at the GeoJSON spec, the duplicate vertex seems to be part of the specification, so I would say in your convert function you skip the last element when converting to points to use for an s2.Loop

o A linear ring is a closed LineString with four or more positions.

o The first and last positions are equivalent, and they MUST contain identical values; their representation SHOULD also be identical.

Do you know what version of the GeoJSON these loops are coming from? Older versions were not required to output points in right-hand rule order. The points you list are in clockwise order, and s2.Loop assumes counter-clockwise. So what you are getting is actually 2 loops that encompass 99.99% of the surface of the earth instead of two small loops.

If you print the RectBound for the loops, you will see that the loops as created end up with a full Rect:

bound:{Lat:{Lo:-1.5707963267948966 Hi:1.5707963267948966} Lng:{Lo:-3.141592653589793 Hi:3.141592653589793}} Lat: -pi/2 (-90) - +pi (90) , Lng: -pi (-180) - pi (180)

If the points are reversed to be counter-clockwise order, the bounds now are fairly tight around the shape: e.g.,

bound:{Lat:{Lo:-1.3066650143192748 Hi:-1.294956760479414} Lng:{Lo:0.6941321285173062 Hi:0.7015497792334233}} Lat: -74.8663906 - -74.195557 Lng: 39.770841388572755 - 40.1958414685095135

You will probably want to call Normalize() on the loops which will invert them if they cover more than half the sphere.

l1.Normalize() l2.Normalize() l1.Intersects(l2) == false now

On Wed, May 22, 2019 at 11:02 AM Yanhe Chen notifications@github.com wrote:

@rsned https://github.com/rsned Thanks for answering. Is there an efficient way to make a valid loop then?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/golang/geo/issues/51?email_source=notifications&email_token=ACBHGBDJ6PEKA5F4J3LHDJTPWWDD3A5CNFSM4HOPUNO2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODV72UBA#issuecomment-494905860, or mute the thread https://github.com/notifications/unsubscribe-auth/ACBHGBGMCFE766RLDE46PHTPWWDD3ANCNFSM4HOPUNOQ .

outsidebedisdoor commented 5 years ago

@rsned Thanks for the clarification. However, if the point data is p1: [ [ [-74.93335056304832,39.283340980346445], [-74.93335056305031,40.30000217770297], [-75.91667795181374,40.30000217770097], [-75.91667795181175,40.28334098034444], [-74.93335056304832,39.283340980346445] ] ] p2: [ [ [-78.4608357441789,32.7877842013595], [-79.2830582223977,33.000284223648], [-79.4513915952751,32.7700063882363], [-78.5544468700461,32.4825063507657], [-78.4608357441789,32.7877842013595] ] ]

On the map, those two are not intersects. I removed the last one and it has been validate. Then normalize them. However, they are still considering intersecting.

outsidebedisdoor commented 5 years ago

@rsned The rectbound is : p1: [Lo[30.9799608, 150.2947587], Hi[89.2303567, -153.3647320]] p2: [Lo[61.1105219, -180.0000000], Hi[90.0000000, 180.0000000]]

outsidebedisdoor commented 5 years ago

Some other examples also fails, such as I construct a Polygon using Rectangle Points: minLat: 30, maxLat: 31, minLng: -71, maxLng: -70

func convertRectangleToPolygon(rectangle *models.Rectangle) s2.Loop {
    var Vertices []s2.Point
    Vertices = append(Vertices, s2.PointFromLatLng(s2.LatLng{s1.Angle(*rectangle.MinLat), s1.Angle(*rectangle.MaxLng)}))
    Vertices = append(Vertices, s2.PointFromLatLng(s2.LatLng{s1.Angle(*rectangle.MaxLat), s1.Angle(*rectangle.MaxLng)}))
    Vertices = append(Vertices, s2.PointFromLatLng(s2.LatLng{s1.Angle(*rectangle.MaxLat), s1.Angle(*rectangle.MinLng)}))
    Vertices = append(Vertices, s2.PointFromLatLng(s2.LatLng{s1.Angle(*rectangle.MinLat), s1.Angle(*rectangle.MinLng)}))
    loop := *s2.LoopFromPoints(Vertices)
    return loop
}

and then another points is:

[[[-142.848733550224,65.2780581505263 ,0],
[-142.869406610785,65.2816692577486 ,0],
[-142.906606068389,65.2936136946703 ,0],
[-142.933446911371 ,65.3013914676046 ,0],
 [-142.954876639064 ,65.3061136856359 ,0],
[-142.972775531381 ,65.3086136823949 ,0],
 [-142.984559699584 ,65.3102803472599 ,0],
[-142.992692481551 ,65.3138914566595 ,0],
 [-142.999557204145 ,65.3141692334514 ,0] ,
[-143.002526371116 ,65.3140025659437 ,0],
[-143.002526234783 ,65.1760859060438 ,0] ,
[-142.992948458279 ,65.1780581306049 ,0] ,
[-142.977845128404 ,65.1819470210269 ,0] ,
[-142.964788184438 ,65.1836136902078 ,0] ,
[-142.92530207324 ,65.1852803639764 ,0] ,
[-142.913760407519 ,65.1861136989371 ,0] ,
[-142.907678742146 ,65.1880581448932 ,0] ,
[-142.894071246426 ,65.1930581469995 ,0] ,
[-142.875424585382 ,65.2000025938815 ,0] ,
[-142.858695982345 ,65.2088914853316 ,0] ,
[-142.850253488789 ,65.2158359304428 ,0],
 [-142.847262661947 ,65.2222248196358 ,0],
 [-142.847331559995 ,65.2313914861557 ,0],
 [-142.843687676307 ,65.2369470415047 ,0],
 [-142.822621298931, 65.2497248225114 ,0],
 [-142.812817971782 ,65.2561137128883 ,0],
 [-142.814827977926 ,65.2622248232271 ,0] ,
[-142.817098260534 ,65.2677803785489 ,0] ,
[-142.813569930943, 65.2713914899776, 0] ,
[-142.813032435321 ,65.275835933844 ,0] ,
[-142.817096326717 ,65.2780581560244 ,0] ,
[-142.826978270638 ,65.2775025983354 ,0],
 [-142.848733550224 ,65.2780581505263 ,0]]]

After validate and normalize, those two shows that have intersections.

rsned commented 5 years ago

s1,Angle's are in radians not degrees, so your values are all greater then 2*pi. You need to multiply the s1.Angle by s1.Degree to get them, in the right range.

On Wed, May 22, 2019 at 4:43 PM Yanhe Chen notifications@github.com wrote:

Some other examples also fails, such as I construct a Polygon using Rectangle Points: minLat: 30, maxLat: 31, minLng: -71, maxLng: -70

func convertRectangleToPolygon(rectangle models.Rectangle) s2.Loop { var Vertices []s2.Point Vertices = append(Vertices, s2.PointFromLatLng(s2.LatLng{s1.Angle(rectangle.MinLat), s1.Angle(rectangle.MaxLng)})) Vertices = append(Vertices, s2.PointFromLatLng(s2.LatLng{s1.Angle(rectangle.MaxLat), s1.Angle(rectangle.MaxLng)})) Vertices = append(Vertices, s2.PointFromLatLng(s2.LatLng{s1.Angle(rectangle.MaxLat), s1.Angle(rectangle.MinLng)})) Vertices = append(Vertices, s2.PointFromLatLng(s2.LatLng{s1.Angle(rectangle.MinLat), s1.Angle(rectangle.MinLng)})) loop := s2.LoopFromPoints(Vertices) return loop }

and then another points is:

[[[-142.848733550224,65.2780581505263 ,0], [-142.869406610785,65.2816692577486 ,0], [-142.906606068389,65.2936136946703 ,0], [-142.933446911371 ,65.3013914676046 ,0], [-142.954876639064 ,65.3061136856359 ,0], [-142.972775531381 ,65.3086136823949 ,0], [-142.984559699584 ,65.3102803472599 ,0], [-142.992692481551 ,65.3138914566595 ,0], [-142.999557204145 ,65.3141692334514 ,0] , [-143.002526371116 ,65.3140025659437 ,0], [-143.002526234783 ,65.1760859060438 ,0] , [-142.992948458279 ,65.1780581306049 ,0] , [-142.977845128404 ,65.1819470210269 ,0] , [-142.964788184438 ,65.1836136902078 ,0] , [-142.92530207324 ,65.1852803639764 ,0] , [-142.913760407519 ,65.1861136989371 ,0] , [-142.907678742146 ,65.1880581448932 ,0] , [-142.894071246426 ,65.1930581469995 ,0] , [-142.875424585382 ,65.2000025938815 ,0] , [-142.858695982345 ,65.2088914853316 ,0] , [-142.850253488789 ,65.2158359304428 ,0], [-142.847262661947 ,65.2222248196358 ,0], [-142.847331559995 ,65.2313914861557 ,0], [-142.843687676307 ,65.2369470415047 ,0], [-142.822621298931, 65.2497248225114 ,0], [-142.812817971782 ,65.2561137128883 ,0], [-142.814827977926 ,65.2622248232271 ,0] , [-142.817098260534 ,65.2677803785489 ,0] , [-142.813569930943, 65.2713914899776, 0] , [-142.813032435321 ,65.275835933844 ,0] , [-142.817096326717 ,65.2780581560244 ,0] , [-142.826978270638 ,65.2775025983354 ,0], [-142.848733550224 ,65.2780581505263 ,0]]]

After validate and normalize, those two shows that have intersections.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/golang/geo/issues/51?email_source=notifications&email_token=ACBHGBAD5BDHNQIYV3SAGPDPWXLC7A5CNFSM4HOPUNO2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODWAVEPQ#issuecomment-495014462, or mute the thread https://github.com/notifications/unsubscribe-auth/ACBHGBAVI2TO3DAVU5BKZKTPWXLC7ANCNFSM4HOPUNOQ .

outsidebedisdoor commented 5 years ago

@rsned That really solved my problems. Haven't notice the part that angle is in radians....