You can do this analysis in the "spdep" package. In the relevant neighbor functions, if you use "longlat=TRUE", the function calculates great circle distance and returns kilometers as the distance unit. In the below example you could coerce the resulting distance list object ("dist.list") to a matrix or data.frame however, it is quite efficient calculating summary statistics using lapply.
require(sp)
require(spdep)
# Create SpatialPolygonsDataFrame for 3 squares
poly1 <- Polygons(list(Polygon(matrix(c(15.7,42.3,16.7,42.3,16.7,41.6,15.7,41.6,15.7,42.3),
nrow=5, ncol=2, byrow=TRUE))),"1")
poly2 <- Polygons(list(Polygon(matrix(c(15.7,42.3,16.7,42.3,16.7,41.6,15.7,41.6,15.7,42.3)+0.5,
nrow=5, ncol=2, byrow=TRUE))),"2")
poly3 <- Polygons(list(Polygon(matrix(c(13.8, 45.4, 15.6, 45.4,15.6, 43.7,13.8, 43.7,13.8, 45.4),
nrow=5, ncol=2, byrow=TRUE))),"3")
spolys = SpatialPolygons(list(poly1,poly2,poly3),1:3)
spolys <- SpatialPolygonsDataFrame(spolys, data.frame(ID=sapply(slot(spolys, "polygons"),
function(x) slot(x, "ID"))) )
proj4string(spolys) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
# Centroid coordinates (not used but provided for example)
coords <- coordinates(spolys)
# Create K Nearest Neighbor list
skNN.nb <- knn2nb(knearneigh(coordinates(spolys), longlat=TRUE),
row.names=spolys@data$ID)
# Calculate maximum distance for all linkages
maxDist <- max(unlist(nbdists(skNN.nb, coordinates(spolys), longlat=TRUE)))
# Create spdep distance object
sDist <- dnearneigh(coordinates(spolys), 0, maxDist^2, row.names=spolys@data$ID)
summary(sDist, coordinates(spolys), longlat=TRUE)
# Plot neighbor linkages
plot(spolys, border="grey")
plot(sDist, coordinates(spolys), add=TRUE)
# Create neighbor distance list
( dist.list <- nbdists(sDist, coordinates(spolys), longlat=TRUE) )
# Minimum distance
( dist.min <- lapply(dist.list, FUN=min) )
# Distance coefficient of variation
( dist.cv <- lapply(dist.list, FUN=function(x) { sd(x) / mean(x) } ) )
According to this question, How to calculate distance in meter between QgsPoint objects in PyQgis?, it's easy.
You just need to return a QgsPoint with your function pixel2coord() :
def pixel2coord(x, y):
xp = (pixelWidth * x) + originX + (pixelWidth/2)
yp = (pixelHeight * y) + originY + (pixelHeight /2)
pnt = QgsPoint(xp, yp)
return pnt
Then you will store a list of QgsPoint than you can compare to your vector Point (store in your feat list) like that:
for vpoint in feat :
for rpoint in pntRstList:
vgeometry = QgsGeometry().fromPoint(vpoint)
rgeometry = QgsGeometry().fromPoint(rpoint)
dist=vgeometry.distance(rgeometry)
I'm not sure it's the most efficient way to do this but it should works.
Best Answer
You can define a function to select pixel values and then transform that raster into points in order to use pointDistance() function. Here is a piece of the code I have used: my raster is a dem and my selected point proceed from a shp (you can make it easily from a point like your xy <- c(590,61)). Then I select just get raster cells above 450m. Finally, I use that pointDistance function to add the value of distance to the SpatialPointsDataFrame created by subseting raster.