e-mail from Jim on 2013-1-18:

-----

Regarding calculation of precision estimates for counties, I just resuscitated something I had started working up (with Bob's inspiration) back on the final Friday afternoon of the Nov 2012 BIEN WG meeting at NCEAS. Using GADM2 administrative polygons data (http://www.gadm.org), this calculation works by first reprojecting a given county into whatever UTM zone contains the county centroid, then (just as Bob suggested in his email) computing the distance from the centroid to the most distance point on the county border.

Have a look at the attached CSV of results for US counties, and mull over whether this seems useful. Maybe even sanity-check against your favorite county or two! For those who care, the SQL statement is pasted below my sig, along with a 'utmzone' function that I borrowed from the PostGIS wiki.

Note that I did the whole UTM reprojection bit because the GADM2 data are in a geographic projection, for which planar distance calculations aren't in meaningful units and aren't consistent across latitudes. Although postgis provides some functions that know how to work directly on the WGS 84 ellipsoid and don't require a projection, the exact functions we need for this aren't available.

Some misc things that were in my head as I did this:

1. This precision measure is conservative insofar as it basically describes a circle of radius equal to the precision value, and most counties will only actually occupy a fraction of the circle depending on their shape -- think of long skinny counties, or those which have an offshore island under their jurisdiction. Totally reasonable! But this _is_ just one conception of what "precision" means for county data.

2. In GADM2, counties are the lowest admin level for the US, making things easy. But other countries have additional levels below county/parish, meaning that we'd need to add a step that first unions those smaller polygons into a single polygon for the county/parish.

3. I was a little concerned about the, um, precision of these UTM-based precision estimates. Although UTM is convenient insofar as a local zone can be applied when looking at a given smallish area (like a county), distance is not perfectly preserved in this projection, which suffers from distortion as you approach the edges of a given zone. However, I just did some testing of the UTM distortion using points at roughly opposite east-west ends of Santa Barbara County. SB makes for a decent test case because it's a fairly large county that lies across two UTM zones. Distance between the test points I used is 100833 meters in UTM, and 100782 meters using the more accurate ellipsoid calculation. So I'm pretty confident we could say that our county precision estimates are themselves precise to 1 km (and probably usually much better than that), which I assume is fine?

As an aside, I was curious to see Franklin County, VA, showing up as having one of the worst precision values. Turns out GADM2 counts Franklin County as including the independent city of Franklin located a few hundred km east (http://en.wikipedia.org/wiki/Franklin,_Virginia). What can I say, anything we do is only as good as our inputs.

Cheers,

Jim

-- generate precision estimates for all US counties

-- (requires utmzone function, defined below)

SELECT name_1 AS state, name_2 AS county,

       ROUND(ST_Length(

         ST_Longestline(

           ST_Transform(geom, utmzone(ST_Centroid(geom))),

           ST_Centroid(ST_Transform(geom, utmzone(ST_Centroid(geom))))

         )

       )::numeric/1000, 3) AS error_km

  FROM gadm2

  WHERE name_0='United States'

  ORDER BY state, county;

-- function taken from postgis users wiki

CREATE OR REPLACE FUNCTION utmzone(geometry)

   RETURNS integer AS

 $BODY$

 DECLARE

     geomgeog geometry;

     zone int;

     pref int;

 BEGIN

     geomgeog:= ST_Transform($1,4326);

     IF (ST_Y(geomgeog))>0 THEN

        pref:=32600;

     ELSE

        pref:=32700;

     END IF;

     zone:=floor((ST_X(geomgeog)+180)/6)+1;

     RETURN zone+pref;

 END;

 $BODY$ LANGUAGE 'plpgsql' IMMUTABLE

   COST 100;

-----

e-mail from Jim on 2013-1-18:

-----

See attached, same as before but now with centroid coordinates (after

back-transforming from the local UTM zone to WGS84), along with their

calculated precision in km. Note that these centroids may differ from

what other authorities claim, because of differences in the method

and/or differences in the input polygons.

Revised CSV attached, revised SQL statement below.

Cheers,

Jim

SELECT name_1 AS state, name_2 AS county,

       ST_Y(ST_Transform(ST_Centroid(geom), 4326)) AS latitude,

       ST_X(ST_Transform(ST_Centroid(geom), 4326)) AS longitude,

       ROUND(

         ST_Length(ST_Longestline(geom, ST_Centroid(geom))

       )::numeric/1000, 3) AS error_km

  FROM (

    SELECT name_0, name_1, name_2,

           ST_Transform(geom, utmzone(ST_Centroid(geom))) AS geom

      FROM gadm2

      WHERE name_0='United States') AS gadm2_utm

  ORDER BY state, county

-----