Created
April 26, 2013 00:08
-
-
Save mashbridge/5464228 to your computer and use it in GitHub Desktop.
A basic point-in-polygon test using libkml.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| // Implements a basic point-in-polygon test using libkml. | |
| #include <assert.h> | |
| #include <iostream> | |
| #include <string> | |
| #include "kml/base/vec3.h" | |
| #include "kml/dom.h" | |
| using kmlbase::Vec3; | |
| using kmldom::CoordinatesPtr; | |
| using kmldom::ElementPtr; | |
| using kmldom::LinearRingPtr; | |
| using kmldom::OuterBoundaryIsPtr; | |
| using kmldom::PointPtr; | |
| using kmldom::PolygonPtr; | |
| using std::cout; | |
| using std::endl; | |
| // See http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html | |
| // for details on how the ray casting algorithm works. | |
| bool IsPointInPolygon(const PointPtr& point, const PolygonPtr& polygon) { | |
| // Note: should check that point and polygon are well-formed, have | |
| // coordinates, etc. | |
| const Vec3 pt_vec = point->get_coordinates()->get_coordinates_array_at(0); | |
| double p_lon = pt_vec.get_longitude(); | |
| double p_lat = pt_vec.get_latitude(); | |
| const OuterBoundaryIsPtr& outer = polygon->get_outerboundaryis(); | |
| const LinearRingPtr& ring = outer->get_linearring(); | |
| const CoordinatesPtr& ring_coords = ring->get_coordinates(); | |
| bool is_contained = false; | |
| size_t n = ring_coords->get_coordinates_array_size(); | |
| for (size_t i = 0, j = n - 1; i < n; j = i++) { | |
| Vec3 vi = ring_coords->get_coordinates_array_at(i); | |
| double i_lat = vi.get_latitude(); | |
| double i_lon = vi.get_longitude(); | |
| Vec3 vj = ring_coords->get_coordinates_array_at(j); | |
| double j_lat = vj.get_latitude(); | |
| double j_lon = vj.get_longitude(); | |
| if ( ((i_lat > p_lat) != (j_lat > p_lat)) && | |
| (p_lon < (j_lon - i_lon) * (p_lat - i_lat) / (j_lat - i_lat) + i_lat)) { | |
| is_contained = !is_contained; | |
| } | |
| } | |
| return is_contained; | |
| } | |
| int main(int argc, char** argv) { | |
| // Square polygon bounded +/- 1°. | |
| const std::string polygon_str( | |
| "<Polygon>" | |
| "<outerBoundaryIs>" | |
| "<LinearRing>" | |
| "<coordinates>" | |
| "-1, -1, 0" | |
| "1, -1, 0" | |
| "1, 1, 0" | |
| "-1, 1, 0" | |
| "-1, -1, 0" | |
| "</coordinates>" | |
| "</LinearRing>" | |
| "</outerBoundaryIs>" | |
| "</Polygon>"); | |
| ElementPtr polygon_root = kmldom::Parse(polygon_str, NULL); | |
| const PolygonPtr polygon = kmldom::AsPolygon(polygon_root); | |
| const std::string inside( | |
| "<Point>" | |
| "<coordinates>0, 0, 0</coordinates>" | |
| "</Point>"); | |
| ElementPtr point_inside = kmldom::Parse(inside, NULL); | |
| assert(IsPointInPolygon(kmldom::AsPoint(point_inside), polygon)); | |
| const std::string outside( | |
| "<Point>" | |
| "<coordinates>2, 0, 0</coordinates>" | |
| "</Point>"); | |
| ElementPtr point_outside = kmldom::Parse(outside, NULL); | |
| assert(!IsPointInPolygon(kmldom::AsPoint(point_outside), polygon)); | |
| // Points on the edge are contained. | |
| const std::string edge( | |
| "<Point>" | |
| "<coordinates>-1, -1, 0</coordinates>" | |
| "</Point>"); | |
| ElementPtr point_edge = kmldom::Parse(edge, NULL); | |
| assert(IsPointInPolygon(kmldom::AsPoint(point_edge), polygon)); | |
| } |
also, in your example, I suspect the coordinate strings are not parsable properly. I think it is safer to format them like this:
const std::string polygon_str(""
""
""
""
"-1,-1,0 1,-1,0 1,1,0 -1,1,0 -1,-1,0"
""
""
""
"");
@geoffhendrey - There was a small typo and I fixed it. Here is the source code if you are interested:
https://github.com/gumdal/libkml-pointinpolygon
The master branch is untouched, look into the other branch for changes.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I think your function for ray tracing contains an error:
Here is the correction:
if (((i_lat > p_lat) != (j_lat > p_lat))
&& p_lon < (j_lon - i_lon) * (p_lat - i_lat) / (j_lat - i_lat) + i_lon){
Or equivalently
if (((i_lat > p_lat) != (j_lat > p_lat))
&& (p_lon < (i_lon - j_lon) * (p_lat - j_lat) / (i_lat - j_lat)) + j_lon) {
(The error is that you were adding latitude instead of longitide.. You are casting a horizontal ray, therefore you should be offsetting longitude not latitude. If you test on a KML file of USA state boundaries and test some locations you will find the function fails properly detect point-in-polygon. Either of the changes above will correct it).