-
Notifications
You must be signed in to change notification settings - Fork 214
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Incorrect results for geometries or inputs in geographic coordinate system at the edges longitude-latitude range #1279
Comments
@rorviksam thanks for the report. It seems that the problem is the redundant point that represents the pole (note that both points (163 -90) and (75 -90) represent the same point i.e. the south Pole). Is it possible to create a case with a single pole point? (i.e. in your example only one of (163 -90), (75 -90) should appear in the definition of the polygon). In any case this should be fixed. |
@vissarion We can't remove any of the points (163 -90), (75 -90) as it would change the shape of the polygon and will not represent the same polygon. |
@rorviksam But they are the same point i.e. the South pole, actually all points |
For example the following code will return namespace bg = boost::geometry;
namespace bgs = boost::geometry::strategy;
using std::string_literals::operator""s;
using Point = bg::model::point<double, 2, bg::cs::geographic<bg::degree>>;
using Polygon = bg::model::polygon<Point, true>;
Point pt_in_melbourne;
Polygon melbourne, melbourne2;
bg::read_wkt(
"POLYGON ((75 -90, 75 -6, 78 -2, 92 -2, 107 -12, 114.5 -12, 115.14222222222223 "
"-14.136944444444444, 126.05888888888889 -23.396944444444443, 131.84 -21.2025, "
"136.32888888888888 -21.499722222222225, 138.39 -26.225277777777777, 146.53333333333333 "
"-29, 151.970694083691 -33.4304025266695, 152.93055555555554 -35.31638888888889, "
"150.32055555555556 -38.18861111111111, 150 -45, 163 -45, 163 -90, 75 -90))"s, melbourne);
bg::read_wkt(
"POLYGON ((75 -90, 75 -6, 78 -2, 92 -2, 107 -12, 114.5 -12, 115.14222222222223 "
"-14.136944444444444, 126.05888888888889 -23.396944444444443, 131.84 -21.2025, "
"136.32888888888888 -21.499722222222225, 138.39 -26.225277777777777, 146.53333333333333 "
"-29, 151.970694083691 -33.4304025266695, 152.93055555555554 -35.31638888888889, "
"150.32055555555556 -38.18861111111111, 150 -45, 163 -45, 75 -90))"s, melbourne2);
bg::read_wkt("POINT(119 -46.0250000000001)"s, pt_in_melbourne);
bool const covers = bg::covered_by(pt_in_melbourne, melbourne,
bgs::within::geographic_winding<bgs::vincenty>());
bool const covers2 = bg::covered_by(pt_in_melbourne, melbourne2,
bgs::within::geographic_winding<bgs::vincenty>());
std::cout << "eq=" << bg::equals(melbourne2, melbourne) << std::endl; |
@vissarion Though the point is in both polygons, they don't describe the same shape. You can easily visualize this using a tool like dbeaver running this SQL statements. select 'p1', st_geometryfromtext('POINT(119 -46.0250000000001)')
union all
select 'poly1', st_geometryfromtext(concat('POLYGON ((75 -90, 75 -6, 78 -2, 92 -2, 107 -12, 114.5 -12, ',
'115.14222222222223 -14.136944444444444, 126.05888888888889 -23.396944444444443, 131.84 -21.2025, ',
'136.32888888888888 -21.499722222222225, 138.39 -26.225277777777777, 146.53333333333333 -29, ',
'151.970694083691 -33.4304025266695, 152.93055555555554 -35.31638888888889, ',
'150.32055555555556 -38.18861111111111, 150 -45, 163 -45, 163 -90, 75 -90))'))
union all
select 'poly2', st_geometryfromtext(concat('POLYGON ((75 -90, 75 -6, 78 -2, 92 -2, 107 -12, 114.5 -12, ',
'115.14222222222223 -14.136944444444444, 126.05888888888889 -23.396944444444443, 131.84 -21.2025, ',
'136.32888888888888 -21.499722222222225, 138.39 -26.225277777777777, 146.53333333333333 -29, ',
'151.970694083691 -33.4304025266695, 152.93055555555554 -35.31638888888889, 150.32055555555556 -38.18861111111111, ',
'150 -45, 163 -45, 75 -90))')); |
What you are visualizing is in cartesian coordinates (i.e. you are using the lon,lat geographic coordinates directly as cartesian coordinates on the Web Mercator Projection---epsg:4326, which is in general wrong because you have to project the geographic polygon to get the "projected" cartesian coordinates) that is why you are creating two different polygons. Those polygons are different indeed but they are not related to the polygons of the issue filed which use geographic coordinates (i.e. they live on a globe). |
@vissarion Accurate thanks for drawing my attention to this. |
@vissarion @barendgehrels This is a continuation of issue #1161 which was only partially fixed.
The code below fails if the geometry contains latitude -90 or 90. Tested on boost 1.85. It will work if you replace -90 by -89.9 or 90 by 89.9.
Geometries with longitude 180 or -180 now seem to work for the few I have tested.
Going deeper into this shows that on boost version 1.82 replacing latitude 90 or -90 by 90 - 1e10 or -90+1e10 provided good results but now you need to drop further in precision to workaround this issue by replacing it with 90 - 1e4 or -90+1e4 on boost version >= 1.83.
The text was updated successfully, but these errors were encountered: