1

I have a shapefile with 50+ different polygonal shapes (representing 50+ different regions) and 10,000+ data points that are supposed to be present in one of the regions. The thing is, the 10,000+ points are already coded with a region they are supposed to be in, and I want to figure out how far they are from this coded region in geo-spatial distance.

My current approach (code below), which involves converting shapefiles to owin objects from the sp library and using distfun gets me distances in lat,long euclidean space. But I would like to get geo-spatial distances (eventually to convert to km). Where should I go next?

#basically cribbed from http://cran.r-project.org/web/packages/spatstat/vignettes/shapefiles.pdf (page 9)
shp <- readShapeSpatial("myShapeFile.shp", proj4string=CRS("+proj=longlat +datum=WGS84"))
regions <- lapply(slot(shp, "polygons"), function(x) SpatialPolygons(list(x)))
windows <- lapply(regions, as.owin)

# need to convert this to geo distance
distance_from_region <- function(regionData, regionName) {
    w <- windows[[regionName]]
    regionData$dists <- distfun(w)(regionData$lat, regionData$long)
    regionData
}
prabhasp
  • 528
  • 3
  • 18
  • Try spDistsN1 in sp, with longlat = TRUE – mdsumner Dec 20 '11 at 20:11
  • @mdsummer: spDistsN1 calculates distances between two points. I need distance from polygon to point. It is not immediately obvious to me how to do the conversion... – prabhasp Dec 20 '11 at 21:47
  • from which coordinates of the polygon? The centroid or some chosen boundary point (the nearest one?), or all of the boundary points? Your question is not clear on this. spDists/spDistsN1 provides the great circle distances you need between sets of points (not between just two points) without choosing a map projection, so it's ultimately what I would use once the right coordinates or summary is extracted - can you clarify the question on this? – mdsumner Dec 20 '11 at 22:00
  • 1
    Sorry, I assumed distance from polygon implied minimal distance = distance from point to nearest point on polygon. I guess I can do distance from a given point to all boundary points in the polygon and then take the minimum... seems a bit inefficient, but I guess workable. Will look for better solution in the meantime. – prabhasp Dec 21 '11 at 06:08

1 Answers1

3

I'd project the data to a euclidean (or near euclidean) coordinate system - unless you are spanning a large chunk of the globe then this is feasible. Use spTransform from maptools or sp or rgdal (I forget which) and convert to a UTM zone near your data.

You also might do better with package rgeos and the gDistance function:

 gDistance by default returns the cartesian minimum distance
 between the two geometries in the units of the current projection.

If your data is over a large chunk of globe then... tricky... 42...

Barry

Spacedman
  • 92,590
  • 12
  • 140
  • 224