I have a dataframe with 3 columns , which has my custom lat / long with a failure rate associated with each Lat / long


        Latitude        Longitude Failure.rate
         51.3           -0.108             0
         51.3           -0.109             0
         51.3           -0.121             0
         51.3            0.0632            1
         51.3            0.0312            0
         51.3           -0.126             0

I have a Polygon shape file which consist of grids over a region. I want to combine both the shape file and point file in R to be able to identify which points + failure rate is associated with which grid and create a shp file as a result to be plotted in any visualisation tool

I tried both over and poly function in R, But i get an empty file filled with NA's . Please find the polygon shape file here

sample aggregated csv ( point files ) I believe it is to do with the different projections used in both point and polygon file , Can someone help me to fix the same. Please find the script below


#Read the points coordinates from the CSV file into a dataframe and convert it to SpatialPointsDataFrame
ID <- read.csv("df_custom_long_lat.csv") 

# Coerce into SpatialPointsDataframe from the data-frame 
coordinates(ID) <- c("Longitude", "Latitude")

#read in the spatial polygon shape file from current directory 
#LGA <- readOGR(".", "grid_shape_file") 

#have a sneak peak of what is loaded in the data , Uses a different geocode since the co-ordinates are not interpretable in relation to points df 
test_LGA <-as.data.frame(LGA@data)

#Assign Similar Projection's -- Not sure which projections are being assigned 
proj4string(ID) <- proj4string(LGA)
#proj4string(ID) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
#proj4string(ID) # to check projections

#Find the point in the polygon 
#over(points , polygon )  # Over = at the spatial locations of object x,
# retrieves the indexes or attributes from spatial object y

inside.LGA <- !is.na(over(ID, as(LGA, "SpatialPolygons"))) 
ID$LGA <- over(ID, LGA)$LGA

write.csv(ID, "ID.csv", row.names=FALSE)
writeOGR(ID, ".", "ID", driver="ESRI Shapefile")

# Using alternate function to over (sp)
#Spatial Eco Package point in polygon , The crs layer should be similar for both points an polygons

test_a <- point.in.poly(ID,LGA,sp = TRUE , duplicate = TRUE)
test_a1 <-as.data.frame(test_a@data)
head(test_a1) # all NA's, failed to assosiate polygons to point dataframes
summary(test_a1) # Index is null 
  • Look at the proj4string of ID and try plotting the two features on top of each other. You cannot simply assign a projection, the data must be transformed into a naw projection space. You will likely need to assin a geographic string to your points and then use spTransform to project it to match the ID grid. May 30 '19 at 14:32
  • @JeffreyEvans , I tried assigning a different projections using the below however how to be sure which is the correct projection #proj4string(ID) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
    – EricA
    May 30 '19 at 14:38
  • Can you recommend a projection to be used especially for the shape file, Apologies i have no experience in spatial geometries
    – EricA
    May 30 '19 at 14:39
  • If your coordinates are, in fact, lat/long then the above string is correct. May 30 '19 at 15:00
  • @JeffreyEvans , I still get output from over function as NA's
    – EricA
    May 30 '19 at 15:10

You are having an issue with only a small subset of your points actually intersecting your polygon grid as well as the need for a projection transformation from geographic to Mercator.

It is difficult to evaluate the results because all the attributes associated with the grid, that occur outside the intersection, will be nodata (NA). In just looking at the data it seems that there are no results being returned but, this is in fact not the case. Out of n=136,210 point locations only n=39,157 occur within the London polygon grid.

Here is some code for evaluating the intersection, using a plot, and then subsetting the points to only those that intersect. I changed the input data and attribute names according to the data that you posted on GoogleDrive and not what is in your post.

First, lets import libraries, add data and transform the geographic point data to a Mercator projection, to match the LGA London grid.



LGA <- rgdal::readOGR(".", "london") 

ID <- read.csv("calculate_data_long_lat.csv") 
  sp::coordinates(ID) <- c("LocationLongitude", "LocationLatitude")
  sp::proj4string(ID) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

ID <- sp::spTransform(ID, proj4string(LGA))

Now we can plot first 100 points and add London grid to plot (takes some time). Notice that the grid does not show up.

plot(ID[1:100,], pch=19, col="red")
  plot(LGA, add=TRUE)

We can do this much more efficiently by using extents. This represents the point and polygon vector data as simple "boxes" showing the maximum extents of the sp objects. After seeing the limit of the London grid within the points extent, we can create a plot that includes the points.

# Plot just extents
plot( as(raster::extent(ID), "SpatialPolygons"))
  plot( as(raster::extent(LGA), "SpatialPolygons"), border="red", add=TRUE )

# Plot extents with points
plot(as(raster::extent(ID), "SpatialPolygons"))
  points(ID, pch=19, col="red", cex=0.50)
  plot( as(raster::extent(LGA), "SpatialPolygons"), add=TRUE )

Now that we have validated the intersection and have an idea of the expected overlap we can intersect points with polygon grid. We can then check to see how many observations fall within the London polygon grid.

ID_LGA <- spatialEco::point.in.poly(ID, LGA, sp = TRUE, duplicate = TRUE)
  cat("Total number of point obs:", nrow(ID), "\n")
  cat("Number of point obs that fall within London grid:",
       length(which(!is.na(ID_LGA$Index))), "\n")

Now, lets subset the intersecting points to remove NA values and them and look at first 10 rows

ID_LGA <- ID_LGA[which(!is.na(ID_LGA$Index)),]

Low and behold, we now have attributes from the London polygon grid. For a gut check, lets plot the results just to verify the overlap of the subset.

  points(ID_LGA, pch=19, col="red", cex=0.50)
  • @Jeffery , Thank you it works smooth and i can understand more of the underlying concepts and able to tie all of them nicely. Thank you for pointers and explanation. Would you be able to refer to decent documentation for people starting in geospatial analysis
    – EricA
    May 30 '19 at 20:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.