1

I am using gDistance function from rgeos package to calculate the nearest distance from points to a polygon. Originally I used for loop to calculate the distance for each row in my points data set, however, the speed is extremely slow, especially when I have a large data set. So I tried to use 'apply'. My codes look like:

Invasive_dist <- apply(sp_invasive_pixels, 1, function(x) gDistance(x, sp_tumor))

whereas sp_invasive_pixels is my SpatialPoint object, and sp_tumor is the SpatialPolygons object. My code does not work, the error message is:

Error in apply(sp_invasive_pixels, 1, function(x) gDistance(x, sp_tumor)) : dim(X) must have a positive length

What do I need to change?

test program:

test_poly <- data.frame(matrix(nrow = 4, ncol = 2))
test_pts <- data.frame(matrix(nrow = 4, ncol = 2))

test_poly$X1 <- c(1, 1, 2, 2)
test_poly$X2 <- c(1, 2, 1, 2)

test_pts$X1 <- c(4, 5, 6, 7)
test_pts$X2 <- c(5, 6, 7, 8)


sp_pts <- SpatialPoints(test_pts)

sp_poly <- SpatialPolygons(list(Polygons(list(Polygon(test_poly)),1)))
  • Did you look at the byid option? Can you make a simple reproducible example using demo data that we can all run that illustrates your problem? – Spacedman Aug 9 at 18:29
  • @Spacedman yes I can, but how to upload the R data? – HAOYANG Mi Aug 9 at 18:35
  • You don't have to upload your data - if its a polygon and some points you should be able to make an illustration using sample data from R packages with maybe test points created from data frames made in a few lines of code - we might only need one polygon and three points to show you how its done. – Spacedman Aug 9 at 18:44
  • @Spacedman Sure, I will edit the question – HAOYANG Mi Aug 9 at 18:46
  • I have been using columbus <- raster::shapefile(system.file("etc/shapes/columbus.shp",package="spdep")) (needs raster and spdep packages installed) to get a set of polygons, extract the first one mydata = columbus[1,] if you only have one in your data. – Spacedman Aug 9 at 18:47
3

Using your test data, this:

> gDistance(sp_poly, sp_pts)
[1] 3.605551

returns the distance from the polygon to the nearest point, but if you add byid=TRUE:

> gDistance(sp_poly, sp_pts, byid=TRUE)
         1
1 3.605551
2 5.000000
3 6.403124
4 7.810250

you get a matrix of distances - each row is one of your points. If you have more than one polygon then you get multiple columns.

You can stop reading now if you want...

Let's look at why apply failed for you. First can I reproduce that with your test data:

> apply(sp_pts, 1, function(x){x})
Error in apply(sp_pts, 1, function(x) { : 
  dim(X) must have a positive length

yes, and I've not needed gDistance there so its something more fundamental with the apply and its not gDistance's fault. I suspectXhere is your points, so what'sdim(sp_pts)`:

> dim(sp_pts)
NULL
> 

apply works on data frames, matrices and arrays, but your object has no dimensions!

Perhaps if we make a spatial points data frame - which is a data frame with some geometry it will work:

> sp_ptsdf = SpatialPointsDataFrame(sp_pts, data.frame(name=LETTERS[1:4]))
> dim(sp_ptsdf)
[1] 4 1
> apply(sp_ptsdf, 1, function(x){x})
Error in as.vector(data) : 
  no method for coercing this S4 class to a vector

no, in this case we get a different error and try another approach. Instead of apply we use sapply to loop over the indexes of your points (in your test it is 1 to 4) and extract each point inside the function:

> sapply(seq_along(sp_pts), function(i){gDistance(sp_pts[i], sp_poly)})
[1] 3.605551 5.000000 6.403124 7.810250

Finally something that works! But note this is the same answer as using byid=TRUE but it will be much slower. I only add this to explain where your error is coming from and to show a general technique which might be useful if you want to apply some other function by rows of a spatial points object.

  • This is a really detailed explanation! thank you so much – HAOYANG Mi Aug 9 at 20:09
  • Also wondering whether use 'byid = TRUE ' is the fastest way, since the point set includes over 500k points, so that the speed is really important. – HAOYANG Mi Aug 9 at 20:44
  • Try it. Also, try using the sf package for your point and polygon objects, and st_distance. – Spacedman Aug 10 at 16:54

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.