Point in polygon calculation using vector geometric methods with application to geospatial data. (2024)

Our proposed algorithm is based on the ray casting method. Suppose we wantto determine whether a point $P(x,y)$ lies within a closed polygon$R$ with vertices $R_{0},R_{1},\ldots,R_{n}$, where $R_{0}=R_{n}.$In the ray casting method, a ray $l$ is drawn moving out from thepoint $P$ in a specific direction. The point $P$ lies inside the polygon$R$ if $l$ intersects the edges of $R$ an odd number of times.In vector geometry, one way of writing out the equation of a lineis the normal form. Let $\boldsymbol{n}$ be a normal vector to theline. Note that a line in the 2D plane has two normals going in oppositedirections that can be denoted $\boldsymbol{n}_{1}$ and $\boldsymbol{n}_{2}$with $\boldsymbol{n}_{1}=-\boldsymbol{n}_{2},$ so $\boldsymbol{n}$can be either $\boldsymbol{n}_{1}$ or $\boldsymbol{n}_{2}.$ Assuming$\boldsymbol{a}$ is a known point on the line and $\boldsymbol{r}$is any point on the line, the equation of the line can be writtenas

Every point on the line satisfies equation 1.Points that do not lie on the line do not satisfy the equation, hence$\left(\boldsymbol{r}-\boldsymbol{a}\right)\cdot\boldsymbol{n}\neq 0$with points on opposite sides of the line having different signs.This is illustrated in Figure 1. The point $C$ is onone side of the line while the point $D$ is on the opposite side.Using the concept of perpendicular distance, the position vectorsof the points $C$ and $D$ can be written as $\boldsymbol{c}=\boldsymbol{e}+\lambda\boldsymbol{\hat{n}}_{1}$and $\boldsymbol{d}=\boldsymbol{f}+\mu\boldsymbol{\hat{n}}_{2},$where $\hat{\boldsymbol{n}}$ represents the unit vector in the directionof $\boldsymbol{n},$ and $\lambda$ and $\mu$ are positive constants.Let $\boldsymbol{n}=\boldsymbol{n}_{1}.$ Then $\boldsymbol{c}=\boldsymbol{e}+\lambda\boldsymbol{\hat{n}}$and $\boldsymbol{d}=\boldsymbol{f}-\mu\boldsymbol{\hat{n}}.$ Substituting$\boldsymbol{c}$ and $\boldsymbol{d}$ into the equation of the linegives

Now $\left(\boldsymbol{e}-\boldsymbol{a}\right)\cdot\boldsymbol{n}=0$,since the points $A$ and $E$ lie on the line. Similarly

In our work, we choose to shoot a ray from the point $P$ with directionvector $\boldsymbol{d}=(1,0).$ This is a ray shot horizontally andto the right of the point $P.$ Then the normal vectors to the lineare $\boldsymbol{n}=(0,\pm 1).$ We choose the vector $\boldsymbol{n}=(0,1).$Since the source of the ray is $P,$ the vector $\boldsymbol{p}$is a point on the line. This gives the equation of the line $l:\left(\boldsymbol{r}-\boldsymbol{p}\right)\cdot\boldsymbol{n}=0.$Thus the line is the locus of points $\boldsymbol{r}\left(x,y\right)$such that $\left(\boldsymbol{r}-\boldsymbol{p}\right)\cdot\boldsymbol{n}=0.$This equation reduces to

For the first constraint which is looking for ray crossings, if theray $l$ crosses the polygon $R$ at the edge connecting the vertices$R_{i}$ and $R_{i+1},$ then the points $R_{i}$ and $R_{i+1}$ willbe on opposite sides of the ray $l,$ hence will have opposite signs.Let $f\left(\boldsymbol{r}\right)=\left(\boldsymbol{r}-\boldsymbol{p}\right)\cdot\boldsymbol{n}=r_{y}-p_{y}.$Hence $f\left(\boldsymbol{r}_{i}\right)\cdot f\left(\boldsymbol{r}_{i+1}\right)\leq 0$if the ray crosses the edge $R_{i},R_{i+1}$, otherwise $f\left(\boldsymbol{r}_{i}\right)\cdot f\left(\boldsymbol{r}_{i+1}\right)>0.$This test will however test for all edge crossings of the ray to theleft and right of $P$.

The second constraint of the ray being cast towards the right fromthe point $P$ will be enforced by testing for the vertices with $x-$componentsgreater or equal to the $x-$component of the point $P;$ $r_{x_{i}}\geq p_{x}.$We tested this contraint in two ways:

1. 1.

   Assume the points $R_{i}$ and $R_{i+1}$ are the endpoints of theedge where the ray cast from $P$ intersects the polygon $R.$ Theconstraint is satisfied if $P$ is on the left of the line $R_{i}R_{i+1}.$To test which side of the line the point $P$ lies, define the directionvector of the line $\boldsymbol{d}=\boldsymbol{r}_{i+1}-\boldsymbol{r}_{i}=(d_{x},d_{y}),$then $\boldsymbol{n}=\left(-d_{y},d_{x}\right)$ or $\boldsymbol{n}=\left(d_{y},-d_{x}\right).$Since we want to test which side of the line $P$ lies on, we choose$\boldsymbol{n}$ to be the vector that forms an acute angle withthe vector $\left(1,0\right).$ We obtain the equation of the lineto be $l_{2}:\left(\boldsymbol{r}-\boldsymbol{r}_{i}\right)\cdot\boldsymbol{n}=0.$Note that the vector $\boldsymbol{r}_{i}$ can be replaced by thevector $\boldsymbol{r}_{i+1}$ in the equation of the line $l_{2}.$Then the point lies to the left of the edge if $l_{2}:\left(\boldsymbol{p}-\boldsymbol{r}_{i}\right)\cdot\boldsymbol{n}\leq 0.$

2. 2.

   The vector form of the equation of the line is given as $\boldsymbol{r}=\boldsymbol{r}_{i}+\lambda\boldsymbol{d}.$This means that

   	$\displaystyle(x,y)$	$\displaystyle=(x_{i},y_{i})+\lambda(d_{x},d_{y})$		(2)
   	$\displaystyle x$	$\displaystyle=x_{i}+\lambda d_{x}$		(3)
   	$\displaystyle y$	$\displaystyle=y_{i}+\lambda d_{y}.$		(4)

   Since the ray is a horizontal line passing through the point $P$we know that the $y-$component of the point of intersection betweenthe edge and the ray is $p_{y}.$ This gives

   $$\lambda=\dfrac{p_{y}-y_{i}}{d_{y}}.$$

   (5)

   We then substitute $\lambda$ to get $x=x_{i}+\lambda d_{x}.$ Theray is then understood to cross the edge if $x\geq p_{x}.$

This helps us to develop the following algorithm:

1. 1.

   Input: $p_{x},$ $p_{y},$ vertices of the polygon $R$ as a 2-dimensionalarray.

2. 2.

   Find the values $f_{i}=r_{y}^{i}-p_{y}$ where $r_{y}^{i}$ is the$y-$component of $R_{i}.$

3. 3.

   Calculate sgn_chg${}_{i}=f_{i}\cdot f_{i+1}.$

4. 4.

   Find $i$ where sgn_chg${}_{i}\leq 0.$ This means there is a ray crossingbetween the vertices $R_{i}$ and $R_{i+1}.$

   See Also

   Optimizing Performance in Numba: Advanced Techniques for Parallelization - GeeksforGeeksBrandon Rohrer on LinkedIn: Numba rule of thumb #5: Use @njit rather than @jit. This tip is already…Brandon Rohrer on LinkedIn: Numba rule of thumb #6: Call Numba-jitted functions once before kicking…Python Tokenizer.sequences_to_matrix Examples, keras.preprocessing.text.Tokenizer.sequences_to_matrix Python Examples

5. 5.

   For each sgn_chg_i found in 4, calculate

   	$\displaystyle\lambda$	$\displaystyle=\dfrac{p_{y}-r_{y}^{i}}{d_{y}}$		(6)
   	$\displaystyle x$	$\displaystyle=r_{x}^{i}+\lambda d_{x}.$		(7)

6. 6.

   Find $c=x\geq p_{x}.$

7. 7.

   If $c$ is odd, then $R$ contains $P.$

Tables 3-8and Figures 6-7show the execution times of the various algorithms for different numberof points from 10 to 17,097,823. It can be seen that generally, shapelyand by extension geopandas are the worst performing methods. In Figure 6(a), the steep slope of the shapely algorithm makes it difficult tosee the performance of the other algorithms. The shapely algorithmis therefore removed in the (b) and (c) parts of the figure to properlyshow the performance of the other algorithms. Figure 6(b) shows the performance of the algorithm for logarithmically spacedpoints from 10–10,000,000 and for 17,097,823 points, while the (c) part shows theperformance of the algorithm for the logarithmically spaced points10, 100 and 1000 points which are not visible in (b). The shapelyexecution time is removed from subsequent plots to allow the executiontimes from the other algorithms to be visualized. It can be seen that,consistently the numba implementation with parallelization performsthe fastest, followed by the numba implementation without parallelization,then the OpenCV implementation, then Algorithm 2, followed by Algorithm1. In Figure 7, the executiontimes for the numba implementations are absent because numba was notsupported on the new python 3.11 as at the time of performing thesimulations. In Tables 7-8,we show the execution times of the algorithms on a linearly spacedscale from 10–10,000 in Table 7 andfrom 1,908,647 to 17,097,823 in Table 8.We perform a linear least squares fit of the data in Table 7to find the equation of the form $ax+b=y$, where $x$ is the numberof points and $y$ is the execution time using the matrix equation

In Equation 8, $\boldsymbol{x}$ is the vector of the number of points thatwere tested and $\boldsymbol{y}$ is the vector of execution times.$\boldsymbol{1}$ is the vector of ones with the same length as $\boldsymbol{x}$and $\boldsymbol{y}.$ The slopes and intercepts of the least squareslines of best fit are shown in Table 9. Figure 8shows the least squares lines and the original data used to calculatethem. Figure 9 shows the least squareslines overlaid on the execution times in Table 8to verify the linearity of the algorithms. It can be seen that thelargest deviations in execution times come from the numba algorithmwith parallelization which may be due to the calls to the multipleprocessor threads. Tables 10 and 11show the absolute and relative errors respectively in the executiontimes projected by the least squares lines for the points in Table8 and the actual execution times inthe table. We can say that it is possible to use the least squaresfit of the execution times in Table 7with between 10 and 10,000 points to predict the execution times fora larger number of points as contained in Table 8.In our case, the largest errors in absolute terms were in the opencv implementation 44.48 seconds error while the numba implementation with parallelization had the highest relative error. However and the least error was in the numba implementation without parallelization which was 11 seconds.