# Calculate longest distance, in KM, between county centroids and their furthest point with PostGIS

Using this shapefile of the US counties at 1:500,000 resolution, I need to know, per county, the length from the centroid of that county to its furthest point away on its polygon. I would like the result in kilometers.

I have uploaded the shapefile to a PostGIS table with:

``````shp2pgsql -W LATIN1 gz_2010_us_050_00_500k.shp public.county_shapes | psql -h localhost -d counties -U jerickson
``````

I did not know which SRID to use so I did not specify one.

Then, modifying this question/answer, I tried the following query, limiting to `name = 'Missoula` (a county in Montana that I am familiar with) for now for testing:

``````SELECT
t.gid AS county_id,
max(ST_Distance(dump.geom, ST_Centroid(t.geom))) AS distance
FROM county_shapes t
JOIN ST_DumpPoints(t.geom) dump ON true
where name = 'Missoula'
GROUP BY 1
ORDER BY 1;
``````

While this seems to be correct (logically), I am not sure how to get the resulting value (`0.9038766878461234`) into kilometers. Using Google Earth to approximate the value, I would expect it to be at least 75 km.

How can I modify my process to get the distance from the centroid to the farthest point on the county polygon for each county in the shapefile?

## 1 Answer

The value you get there is most likely in decimal degrees...and thus quite useless; a degree doesn't represent the same ground distance in different latitudes and getting from there to metric values is tricky and unnecessary.

First off, you should always define a SRID for your geometries. You can, for example, run ogrinfo on the shapefile, the CRS is usually defined in the metadata.
Now, since ST_Distance returned decimal degrees, it's definetely referenced in a geographical CRS. Possibly the direct and most precise way then is to cast your geometries to geography type, i.e.:

``````... max(ST_Distance(dump.geom::geography, ST_Centroid(t.geom)::geography)) AS distance ...
``````

This should return values in meter, based on an unknown datum, but calculated on a spheroid and possibly still quite precise, despite the missing CRS. Which is no excuse for not having tried harder to find one!

Update for completeness:
Use `... max(...) / 1000 ...` to get results in km.

• Thank you, this is extremely helpful. Casting to `geography` type definitely provides a value that matches my estimate. When running `ogrinfo` on that shapefile, I get only `1: gz_2010_us_050_00_500k (Polygon)`. Is there a SRID that is accepted for use with the US as a whole? Or would one need to do it state-by-state, for example? – Jeff Erickson Mar 20 '18 at 21:58
• @JeffErickson nah, you can't just assume one ,) use `ogrinfo <shapefile> gz_2010_us_050_00_500k -so` to access the actual layers metadata. – ThingumaBob Mar 21 '18 at 8:19
• Thank you. The output is `Layer SRS WKT: GEOGCS["GCS_North_American_1983", DATUM["North_American_Datum_1983", SPHEROID["GRS_1980",6378137,298.257222101]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]]`. Does this mean we are using EPSG:4269? – Jeff Erickson Mar 22 '18 at 12:20
• @JeffErickson that's the one! – ThingumaBob Mar 22 '18 at 12:33