3

Very simple situation : a polygon define a geographical area and I want to know whether a point, given by it gps coordinates, lies within that polygon.

I went through many SO questions and have tried various functions and packages like sp, but cannot make out why it fails.

I tried with this very simple function: https://www.rdocumentation.org/packages/SDMTools/versions/1.1-221/topics/pnt.in.poly

install.packages("SDMTools v1.1-221")
library(SDMTools v1.1-221)

## Coordinates of the polygon corners
lat <- c(48.43119, 48.43119, 48.42647, 48.400031, 48.39775, 48.40624, 48.42060, 48.42544, 48.42943 ) 
lon <- c(-71.06970, -71.04180, -71.03889, -71.04944, -71.05991, -71.06764, -71.06223, -71.06987, -71.07004)
pol = cbind(lat=lat,lng=lon)

## Point to be tested
x <- data.frame(lng=-71.05609, lat=48.40909)

## Visualization, this point clearly stands in the middle of the polygon
plot(rbind(pol, x))
polygon(pol,col='#99999990')

## Is that point in the polygon? 
out = pnt.in.poly(x,poly)

## Well, no (pip=0)
print(out)

The example given for this function works with me, but this simple case no... why is that?

Xavier Prudent
  • 1,570
  • 3
  • 25
  • 54

2 Answers2

5

I have not used the method that you are using, but I have one from within sp which works flawlessly on your point and polygon.

I cherry picked your code and left the lat and lon as vectors and the point coordinates as values to suit the functions requirements.

But you could just has easily have made a data frame and used the columns explicitly as lat/lon values.

Here is the gist of it:

require(sp)
## Your polygon
lat <- c(48.43119, 48.43119, 48.42647, 48.400031, 48.39775, 48.40624, 48.42060, 48.42544, 48.42943 ) 
lon <- c(-71.06970, -71.04180, -71.03889, -71.04944, -71.05991, -71.06764, -71.06223, -71.06987, -71.07004)


## Your Point
lng=-71.05609
lt=48.40909

# sp function which tests for points in polygons

point.in.polygon(lt, lng, lat, lon, mode.checked=FALSE)

And here is the output:

[1] 1

The interpretation of this from the documentation:

integer array values are:

  • 0 point is strictly exterior to polygon
  • 1 point is strictly interior to polygon
  • 2 point lies on the relative interior of an edge of polygon
  • 3 point is a vertex of polygon

As your point is a 1 based on this, it should be wholly within the polygon as your map shows! The key to getting good output with these types of data is serving the variables in the right formats.

you could just as easily had a data frame df with df$lat and df$lon as the two polygon variables as well as a test frame test with test$lat and test$lon as a series of points. You would just substitute each of those in the equation as such:

point.in.polygon(df$lat, df$lon, test$lat, test$lon, mode.checked=FALSE)

And it would return a vector of 0's, 1's 2's and 3's

Just be sure you get it in the right format first! Here is a link to the function page:

sconfluentus
  • 4,693
  • 1
  • 21
  • 40
2

I can't see it explicitly stated in the documentation for ?pnt.in.poly, but it appears the ordering of the lng and lat columns matter. You need to swap the column ordering in your pol and it works. 

pol = cbind(lat=lat, lng=lon)
pnt.in.poly(x, pol)
#          lng      lat pip
# 1 -71.05609 48.40909   0

pol = cbind(lng=lon, lat=lat)
pnt.in.poly(x, pol)
#          lng      lat pip
# 1 -71.05609 48.40909   1

In spatial goemetry, lng is often thought of as the x-axis, and lat the y-axis, which you'll see is reversed in your plot()

SymbolixAU
  • 25,502
  • 4
  • 67
  • 139