point_in_polygon - ObjectVision/GeoDMS GitHub Wiki

Geometric functions point_in_polygon

• point_in_polygon(point_data_item, polygon_data_item)

point_in_polygon(point_data_item, polygon_data_item) results for each point of the point_data_item in a relation towards the domain unit of the polygon_data_item in which the point is located.

The resulting values unit is the domain unit of the polygon_data_item.

If a point is not located in any polygon, the function results in a null value (the yellow point in the image). If a point is located in multiple polygons, the function results in one of them, and which one is not specified: it is the first one that the function's spatial index finds, which follows from how that index is built and can differ between versions of GeoDMS. A configuration, and a test, must not depend on which one it is. If you want to choose among overlapping polygons, use the point_in_ranked_polygon function; with a constant rank, that takes the polygon with the lowest index.

The first one found is kept because it is fast: a search that stops at the first polygon containing the point does not have to look at the rest, while the polygon with the lowest index is only known after searching on to the end, which takes longer, also on polygons that do not overlap. A point that lies in only one polygon always gets that polygon.

To rasterize polygon data, use the poly2grid function rather than point_in_polygon.

The function builds a spatial index over the polygons and tests each point against the polygons whose bounding box holds it: O(n log p + n * c), where n = number of points, p = number of polygons, and c = the cost of testing a point against such a polygon. For large-scale rasterization, prefer poly2grid.

The spatial index keeps each polygon at a depth that follows from its size (a loose quadtree), also when the polygon crosses or touches one of the lines along which the index divides its area. So in a dense set of small polygons, such as parcels or grid cells, a point is tested against only the few polygons near it.

A point that is tested against a polygon of many vertices reads only a small part of them. A polygon of at least 32 vertices gets an edge index: its height is cut into horizontal slabs, about one per 8 vertices, each with the edges that reach into it, and a point tests only the edges of its own slab. So c hardly grows with the number of vertices: on one machine, 10 million points over 400 squares take 0.36 s when each square has 10 001 vertices on its boundary, and 0.2 s when it has 5. The edge index does not change the results.

The edge index is built only where it pays off: for the polygons whose bounding box touches the bounding box of all the points, and only when there are at least 4 points per such polygon. So the addresses of one village tested against all neighbourhoods of a country index only the neighbourhoods around the village, and a few points tested against many detailed polygons use no index. The index takes at most 4 edge numbers of 4 bytes per vertex of an indexed polygon, as much as the coordinates of a dpoint polygon, and is released when the function is done. point_in_ranked_polygon and point_in_all_polygons use the same index.

• a data item point_data_item with Point value type
• a data item polygon_data_item with composition type type polygon and Point value type

attribute<district> district_rel (ADomain) := point_in_polygon(Adomain/point, district/geometry);

ADomain, nr of rows = 5

domain District, nr of rows = 7

• point_in_ranked_polygon
• point_in_all_polygons