Sp::Over() for Point in Polygon Analysis

sp::over() for point in polygon analysis

You should not supply a function. You are aggregating the attribute values of your points over the geometry of the polygon, (i.e. the number returned is the mean of the attribute of the points that fall within the polygon). In addition you have your x and y the wrong way round for what you want to do. Should be...

over( pnts , ind_adm , fn = NULL) 

Keep getting NAs when I run Over() function on Points(Lat,Lon) on shapefile polygons

I realized that your point data is an sf object since you have POINT (-82.34323174 29.67058748) as character. Hence, I reconstructed your data first. I assigned a projection here as well.

library(tidyverse)
library(sf)
library(RCurl)

url <- getURL("https://raw.githubusercontent.com/THsTestingGround/SO_readOGR_quest/master/311_Service_Requests__myGNV_.csv")

mydf <- read_csv(url) %>% 
        mutate(Location = gsub(x = Location, pattern = "POINT \\(|\\)", replacement = "")) %>% 
        separate(col = "Location", into = c("lon", "lat"), sep = " ") %>% 
        st_as_sf(coords = c(3,4)) %>% 
        st_set_crs(4326)

I imported your shapefile using sf package since your data (mydf in this demonstration) is an sf object. When I imported the data, I realized that I had LINESTRING, not polygons. I believe this is the reason why over() did not work. Here I created polygons. Specifically, I joined all seven polygons all together.

mypoly <- st_read("cgbound.shp") %>% 
          st_transform(crs = 4326) %>% 
          st_polygonize() %>% 
          st_union()

Let's check how your data points and polygon are like. You surely have data points staying outside of the polygon.

ggplot() +
geom_sf(data = mypoly) +
geom_point(data = mydf, aes(x = Longitude, y = Latitude))

You said, "I need some points which are outside of the polygon to be NA." So I decided to create a new column in mydf using st_intersects(). If a data point stays in the polygon, you see TRUE in the new column, check. Otherwise, you see FALSE.

mutate(mydf,
      check = as.vector(st_intersects(x = mydf, y = mypoly, sparse = FALSE))) -> result

Finally, check how data points are checked.

ggplot() +
geom_sf(data = mypoly) +
geom_point(data = result, aes(x = Longitude, y = Latitude, color = check))

If you wanna use over() mixing with this sf way, you can do the following.

mutate(mydf,
       check = over(as(mydf, "Spatial"), as(mypoly, "Spatial")))

The last thing you wanna do is to subset the data

filter(result, check == TRUE)

THE SIMPLEST WAY

I demonstrated you how things are working with this sf approach. But the following is actually all you need. st_filter() extracts data points staying in mypoly. In this case, data points staying outside are removed. If you do not have to create NAs for these points, this is much easier.

st_filter(x = mydf, y = mypoly, predicate = st_intersects) -> result2

ggplot() +
geom_sf(data = mypoly) +
geom_sf(data = result2)

Point in Polygon method in R for merge regions of a shapefile with over 15'000 observations of a dataset

Here's a somewhat contrived example. The key function is st_intersection(). There may be faster solutions using other methods in sf.

library(sf)
library(tmap)
library(tidyverse)
data("World")
world.centroids <- st_centroid(World)[,1]
world.centroids.expanded <- world.centroids[rep(row.names(world.centroids), 3),]
world.centroids.expanded$ref <- runif(nrow(world.centroids.expanded), 0 ,1 )
intersections <- st_intersection(World[,1], world.centroids.expanded)
intersections <- st_drop_geometry(intersections)
intersections |> group_by(iso_a3) |> summarise(avg = mean(ref))

# A tibble: 164 × 2
   iso_a3   avg
   <fct>  <dbl>
 1 AFG    0.382
 2 AGO    0.457
 3 ALB    0.510
 4 ARE    0.359
 5 ARG    0.627
 6 ARM    0.361
 7 ATA    0.130
 8 ATF    0.241
 9 AUS    0.636
10 AUT    0.304

R keep points within one polygon

I have find something more :-)

pts_in<-pts[!is.na(over(pts,poly)),]

It keep only points in the polygon (source : http://cran.r-project.org)

In R, find the polygon containing a point - Shapefile

You can use point.in.polygon from sp:

library(sp)
library(magrittr)

bairrorio.fort %>% 
  split(.$id) %>% 
  sapply(function(x) point.in.polygon(p[1], p[2], x$long, x$lat) > 0) %>%  
  names(.)[.]

Calculate minimum distance from polygon to spatial point

There are two issues with your current approach.
1.) The coordinates assignment isn't quite correct. It should be long/lat, but you're assigning it as lat/long.
2.) Directly setting the CRS in the current way will not actually alter the points in the necessary way. You need to assign an appropriate long/lat CRS first, and then perform a spTransform action.

#Open data on tin mines and define as spatial points
tin.01 <- read.csv("random/Tin Mine Location_01.csv")
coordinates(tin.01) <- c("Longitude","Latitude")  ## Should be `9:8` if you wanted to stick with indexes, but using the names here is generally lower risk.
proj4string(tin.01) <- CRS(proj4string(Kinta)) ## Setting the initial projection to the same one the polygons are using. You should change this if your original data source uses some other known long/lat projection.

tin.km <- spTransform(tin.01,CRS(EPSG.3375)) ## Creating a transformed set of points for the distance calculation.

#Find distance between district and mines
gDistance(Kinta.km,tin.km,byid=TRUE) ## '0' distance means the mine is inside the district.

          59
1 194384.372
2 223773.999
3      0.000
4  36649.914
5 102944.361
6      0.000
7      0.000
8   6246.066
9      0.000

Find the right polygon contains specific Point

This returns a result:

select AsText(polygon.geom)
from   polygon
where  Contains(polygon.geom, 
           PointFromText('POINT(32.832677661456 39.604395901764)'));

See this fiddle, which has two tables, polygon (with the one value you provided), and location with two records, one of which has the point you provided, and another one with a point outside the given polygon.

The fiddle SQL finds that the first point is within the polygon:

select AsText(polygon.geom), AsText(location.geom)
from   polygon, location
where  Contains(polygon.geom, location.geom);

Polygon and line input for point output (in R)

You can do this with the spatstat package, but you need to convert
from SpatialLines and SpatialPolygons to psp (planar segment
pattern) and owin (observation window -- a polygon in the plane) using
the maptools package.

Since you didn't provide example data I will first construct
SpatialPolygons and SpatialLines objects looking like your graphic:

library(maptools)
#> Loading required package: sp
#> Checking rgeos availability: TRUE
library(spatstat)
#> Loading required package: nlme
#> Loading required package: rpart
#> 
#> spatstat 1.48-0.029       (nickname: 'Alternative Facts') 
#> For an introduction to spatstat, type 'beginner'
A <- owin(c(1,3), c(2,5))
B <- owin(c(4,7), c(1,5))
p <- union.owin(A,B)
p <- as(p, "SpatialPolygons")
l <- psp(c(0.5, 0.5), c(2, 4), c(4, 6), c(2, 4), window = owin(c(0,7), c(0,5)))
l <- as(l, "SpatialLines")
plot(p)
plot(l, col = "green", add = TRUE, lwd = 3)

Now to do what you ask you convert the edges of the polygon to a psp
and find where the lines cross these edges (returned as a planar point
pattern ppp):

l <- as(l, "psp")
p <- as(p, "owin")
e <- edges(p)
x <- crossing.psp(l, e)
x
#> Planar point pattern: 4 points
#> window: rectangle = [1, 6] x [2, 4] units

Adding the points to the plots shows what we found:

plot(p, main = "")
plot(l, col = "green", add = TRUE, lwd = 3)
plot(x, add = TRUE, col = "red", pch = 20, cex = 3)

---

Related Topics