0

Given a bunch of 2D points and a polygon, I want to evaluate which points are on the boundary of the polygon, and which are strictly inside/outside of the polygon.

The 2D points are:

> grp2
          x2         y2
1  -5.233762  1.6213203
2  -1.107843 -7.9349705
3   4.918313  8.9073019
4   7.109651 -3.9571781
5   7.304966 -4.3280168
6   6.080564 -3.5817545
7   8.382685  0.4638735
8   6.812215  6.1610483
9  -4.773094 -3.4260797
10 -3.269638  1.1299852

and the vertices of the polygon are:

> dfC
         px         py
1  7.304966 -4.3280167
2  8.382685  0.4638735
3  6.812215  6.1610483
4  5.854366  7.5499780
5  2.385478  7.0895268
6 -5.233762  1.6213203
7 -4.773094 -3.4260797
8 -1.107843 -7.9349705

The plot of the situation looks like following:enter image description here

Clearly, there are 3 points inside the polygon, 1 point outside and 6 points on the edge (as is also evident from the data points).

Now I am using point.in.polygon to estimate this. According to the documentation of package sp, this should return 'integer array; values are: 0: point is strictly exterior to pol; 1: point is strictly interior to pol; 2: point lies on the relative interior of an edge of pol; 3: point is a vertex of pol.'

But my code is not being able to detect the points which are vertices of the polygon:

> point.in.polygon(grp2$x2,grp2$y2,dfC$px,dfC$py)
 [1] 0 0 0 1 0 1 0 0 0 1

How can I resolve this problem?

Pratik Mullick
  • 55
  • 6
  • I agree. But some of the remaining points are. And these are the points obtained in my code. I need to figure out how can I deal with the 'precision' problem. – Pratik Mullick Aug 10 '22 at 19:00
  • 1
    How precise do you need to be? For example, 7.304966 -4.3280167 on the lower right side is inside the vertex at 7.304966 -4.3280168. One solution is to round appropriately. – John Garland Aug 10 '22 at 19:09
  • @ZheyuanLi ``` > dput(grp2) structure(list(x2 = c(-5.23376158438623, -1.10784274060279, 4.91831264458597, 7.10965098813176, 7.30496606323868, 6.08056389726698, 8.38268484454602, 6.81221520062536, -4.77309438399971, -3.26963831204921), y2 = c(1.6213203035295, -7.93497047852725, 8.90730194281787, -3.95717813633382, -4.32801675051451, -3.58175448607653, 0.463873511180282, 6.16104830056429, -3.42607971746475, 1.12998515367508)), class = "data.frame", row.names = c(NA, -10L )) ``` – Pratik Mullick Aug 10 '22 at 19:14
  • @ZheyuanLi > dput(dfC) structure(list(px = c(7.30496604690398, 8.38268483267204, 6.81221519657874, 5.8543657224258, 2.38547779172217, -5.23376157160271, -4.77309438207759, -1.10784272671463), py = c(-4.32801673416477, 0.463873496893438, 6.16104828907916, 7.5499780328743, 7.0895267949496, 1.6213202900966, -3.42607970231693, -7.9349704726765)), class = "data.frame", row.names = c(NA, -8L)) – Pratik Mullick Aug 10 '22 at 19:14

1 Answers1

0

The points are not equal. For example, grp2$x2[1] == -5.23376158438623 but for fpC$px[6] == -5.23376157160271. These are not equal. As the comments suggest, you will have more luck if you round the values:

grp3 <- round(grp2, 3)
dfC3 <- round(dfC, 3)
point.in.polygon(grp3$x2,grp3$y2,dfC3$px,dfC3$py)
#  [1] 3 3 0 1 3 1 3 3 3 1

Now

grp3[1, ]
#        x2     y2
# 1  -5.234  1.621
fpc3[6, ]
#       px     py
# 6 -5.234  1.621

Changing the number of decimals to 4 or 5 gives the same results as 3. For floating point numbers to be equal, they must match exactly over all 14 decimal places.

dcarlson
  • 10,936
  • 2
  • 15
  • 18