1

I'm trying to calculate pair distances among the points of my dataset. I have point in EPSG:3857.

This is the code I ran:

I import the data, then:

coordinates(datPb) <- c("X","Y")

proj4string(datPb) <- CRS("+init=epsg:3857")

datPb_owin <- as.owin(list(xrange=c(min(datPb$X),max(datPb$X)),yrange=c(min(datPb$Y),max(datPb$Y))))

datPb.ppp <- ppp(x=datPb$X, y=datPb$Y,marks=datPb$Pb, window=datPb_owin)

pairdist(datPb.ppp)

The pairdist function return a matrices of distance but I know (from field measurements and comparison with QGIS on a subsample of points) that these distances are not correct. I though it might be an error in the understanding of the EPSG but when I save back my data as a shapefile and measure distance in QGIS again they are fine.

here is a subset of my data:

quadrat X Y

28 271322.8329 6231012.4650

29 271313.2814 6231006.6400

30 271319.1768 6231013.2850

31 271326.7498 6231017.9710

Improve this question
1
  • When your are defining the projection is it, in fact, the native projection of the coordinates? Assigning a projection is not the same as reprojecting the data using spTransform. If you are defining a projection that is different than what the data is in then this would cause you problems. If you have a shapefile then why are you importing a flat file? Just use rgdal:readOGR to read the shapefile and it will retain the original projection.
    – Jeffrey Evans
    Commented Nov 30, 2016 at 17:06

2 Answers 2

Reset to default
1

It seems like you are over complicating things a bit in using spatstat:pairdist. To return a distance matrix, just use the spDist in sp. This avoids coercing an sp SpatialPointsDataFrame object, which is passed directly to spDist, to a spatstat ppp object. 

library(sp)
pts <- data.frame(ID=c(28,29,30,31),
  x=c(271322.8329, 271313.2814, 271319.1768, 271326.7498),
  y=c(6231012.4650, 6231006.6400, 6231013.2850, 6231017.9710))
coordinates(pts) <- ~x+y
proj4string(pts) <- "+init=epsg:3857"

spDists(pts)

I will add that error needs to be taken into account when evaluating discrepancies between field and in instrument based measurements. In recreation grade GPS units (eg., Garmin) a 10m error is not uncommon. When you calculate the difference between observations [28,30], the stated "field measurement" is 2.49m where the distance calculated from the coordinates sqrt((271322.8329 - 271319.1768)^2 + (6231012.4650 - 6231013.2850)^2) is 3.746928m. Honestly, until you get into survey grade GPS units, a 1.256928m error between two coordinates it pretty good. Although, projection should also be taken into account. Some projections can add quite a bit of distortion. You should check the distance(s) in the native coordinate system to make sure that error is not being introduced by a given projection or datum/spheroid transformation.

Improve this answer
6
  • This gives me the same distance table than with pairdist. Which make me think that it is problem in R in interpreting the projection but I do not understand where is the issue.
    – user41334
    Commented Nov 30, 2016 at 8:32
  • @user41334 I understand that both codes produce the same result. For instance the distance between the pair of #28 - #29 was 11.187572. I am no good at math, but it is probably supported by sqrt((271313.2814 - 271322.8329)^2 + (6231006.6400 - 6231012.4650)^2), too. What in particular makes you think these are not correct? Or number of decimal places ...not enough? Sorry, I am just curious to know what is the issue.
    – Kazuhito
    Commented Nov 30, 2016 at 13:04
  • I know it is wierd but I have measured in the field the distances and for example between pt28 and 30 I should find 2.49. This is the distance that I also find when plotting the points of QGIS. Why is that different, that is the question I was wondering too.
    – user41334
    Commented Nov 30, 2016 at 14:51
  • Not sure on why you are using a Web Mercator but the projection units are in meters and in checking the resulting distances manually, they are correct. What units are your field measurements? What is the error of your GPS? Some consumer grade GPS units can very high vertical error, particularity in mountainous areas.
    – Jeffrey Evans
    Commented Nov 30, 2016 at 16:15
  • you mean that R gives you a distance between 28 and 30 of 2.49m? my field unit are in meters.
    – user41334
    Commented Nov 30, 2016 at 16:39
0

If this might help someone, I found my answer. The issue is that even though you specify the projection of the data in R, the distance matrix you get relies on the unit of the projection and is not transformed in meters (unit). If you wish to get a distance matrix in meters, you should give to R coordinates that are expressing meters because the function it uses to calculate distances is always sqrt(xa-xb)^2+(ya-yb)^2), whichever EPSG projection you specified.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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