It's not a simple problem, it involves some kind of forced iteration over the set of candidate points. This chapter from the workshop shows a similar problem, but not exact (your problem is slightly easier)
http://postgis.net/workshops/postgis-intro/advanced_geometry_construction.html
The nearest neighbor searching chapter from the workshop shows the tools you might use to do an index-assisted approach with some external loop driving the query
http://postgis.net/workshops/postgis-intro/knn.html
If your points have a distinct id and you know a distance tolerance (9999) they will all fall within, a self-join and use of the "DISTINCT ON" filter will get you the answer in one go.
WITH unfiltered AS
(
SELECT t1.id AS id1, t2.id AS id2, ST_Distance(t1.geom, t2.geom) as dist
FROM t t1, t t2 WHERE ST_DWithin(t1.geom, t2.geom, 9999) AND t1.id <> t2.id
ORDER BY t1.id, ST_Distance(t1.geom, t2.geom) ASC
)
SELECT DISTINCT ON (id1) id1, id2, dist FROM unfiltered;
It first gathers the candidates combinations of points, and sorts them by distance. Then the "distinct on" filter strips out just the first member of each candidate group, which conveniently is the closest, thanks to the pre-sorting.
Use gCentroid
from rgeos to calculate centroids.
Use spDists from sp (as mentioned by @Jot eN above) to calculate distance matrices.
I made up a simple example using admin boundaries. Note that the parameters of spDists
change, depending on your map units. The example is for degrees. If you have meters use something like: spDists(centroids)/1000
library(sp)
library(rgdal)
library(rgeos)
# getting sample data of the three main Belgian regions
tmpdir <- tempdir()
download.file("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_map_subunits.zip", "world.zip")
unzip("world.zip", exdir = tmpdir )
shapefile <- paste0(tmpdir,"/ne_10m_admin_0_map_subunits.shp")
world <- readOGR(shapefile, "ne_10m_admin_0_map_subunits")
unlink(tmpdir)
belgium <- world[world$ADM0_A3 == "BEL", ]
centroids <- gCentroid(belgium, byid=TRUE)
spDists(centroids, longlat=TRUE) #distance in km
Best Answer
Since you have ArcInfo, you can use the Near geoprocessing tool to find the distance. It's in ArcToolbox > Analysis > Proximity. Your input features will be Census tract centroids and your near features will be hospitals.
You may want to consider running the analysis simply on the Census tract polygons instead of their centroids. Tracts are large and often very odd-shaped, especially in rural areas. If you use centroids, you can end up with weird-looking results if elongated panhandles are close to a hospital but that tract's centroid is still far.
Moreover, you may also want to consider using a higher spatial resolution Census geography, e.g. block groups or blocks. American Community Survey data are available down to the block group resolution, and would probably produce better looking results.