4

TL;DR:
How do I write a query in SQL (postgres, postGIS, tiger extension) that will return census geographies (blocks, tracts, zcta) within an X mile radius of a lat/lng coordinate?

Background:
I would like to feed a query a lat/lng pair and receive back census geographies (block, blockgroup, tract, county, state).

EDIT Piggybacking off this answer, I have a query that will return the information I need at a single coordinate. How do I extend it to include a radius?
Finding census block for given address using Tiger geocoder

SELECT tabblock_id
FROM tabblock 
WHERE ST_Contains(the_geom, ST_SetSRID(ST_Point(-71.101375, 42.31376), 4269))

This answer here suggests I can use something like this to get the data I need. With the query below and the previous one, I think can mash them up to get the radius information I need.
Return all results within a 30km radius of a specific lat/long point?

SELECT *
FROM your_table
WHERE ST_Distance_Sphere(the_geom, ST_MakePoint(your_lon,your_lat)) <= radius_mi * 1609.34

However, now I'm seeing my tabblock table is empty.

geocoder=# select count(*) from tiger.tabblock;
 count
-------
     0
(1 row)
5
+50

Query structure:

In general, querying a table for true radius lookups is best executed using ST_DWithin, i.e.:

SELECT *
FROM <your_table>
WHERE ST_DWithin(<geom_colum>, <ref_geom>, <distance>);

where <geom_column> refers to the tables geometry column, <ref_geom> the e.g. point geometry in question, and <distance> the radius in units of the CRS (see Notes on CRS). With proper indexes, this should be orders of magnitude faster than other query designs.

If true radius is not of concern, a direct bbox search should again be substantially faster on very large tables, i.e.

SELECT *
FROM <your_table>
WHERE <geom_colum> && ST_Expand(<ref_geom>, <distance>);

This will, of course, find geometries in the rectangular bbox around the <ref_geometry>, and in the case of an expanded point, the square with the <distance> * 2 dimensions.

Note that combining both does not improve execution time in my experience; ST_DWithin already implicitly executes the same bbox filter.


Notes on CRS:

Naturally, of substantial importance is the underlying CRS. ST_DWithin (as well as e.g. ST_Expand) will assume <distance> to be in the units of the data's CRS; if your data is in a geographic projection (e.g. EPSG:4326), units are degrees...and as such useless.

You have two options (pros and cons are from the top of my head):

  • reproject (ST_Transform) your data into a suitable planar projection like State Plane or UTM:

    • pros:

      • working units will be in CRS units, e.g. meter or feet
      • consistent measurements (distance, area) within the limits of the projection
      • generally faster calculations
      • can be reprojected on-the-fly
    • cons:

      • coverage of each system's projections; e.g. with State Plane, you will need to transform each dataset with the respective EPSG, same with UTM zones

      • precision loss due to the limits/distortions of planar projections

  • cast to or directly use the geography data type:

    • pros:

      • working units are in meter
      • working CRS is EPSG:4326; no transformations needed
      • best possible precision of measurements due to sperical/spheroidal algebra
      • can be cast on-the-fly
    • cons:

      • calculations are more demanding; increase of execution time
      • not all functions accept the geography type (e.g. ST_Expand)
      • index scans might not be possible to use


I'm not in the US, the TIGER data and geocoding service are not much of interest to me, and I don´t know the data's details, but I expect the CRS to be EPSG:4326 in order to guarantee nationwide consistency?

If that is the case, I would recommend to use a cast to geography and accept the larger latency in execution, e.g.

SELECT *
FROM <your_table>
WHERE ST_DWithin(<geom_colum>::geography, <ref_geom>::geography, <distance_in_mi> * 1609.34);

or, for completeness, the same with a bbox filter (note the use of _ST_Expand, with a leading _; this 'meta-function' is not exactly encouraged to use, but works on geography type)

SELECT *
FROM <your_table>
WHERE <geom_colum>::geography && _ST_Expand(<ref_geom>::geography, <distance_in_mi> * 1609.34);

Note: you can set the use_spheroid parameter of ST_DWithin to false to increase performance (calculations will be based on a sphere instead of a spheroid; less precision, faster check).

If you intend to run this over and over on a nationwide dataset, consider the extra time and disk space and add a geography column with an own index.


TIGER data (not my expertise):

As to why your tabblock table is empty: did you import the data?

I do hope this actually is you issue; getting geometries from a present table within a radius. About the presence of that table: as far as I got it from a very quick glance at the docs, the TIGER geocoder works on manually pre-loaded data!?

  • Thanks for the detailed explanation. I'm working on setting up the database and getting all the data imported correctly now. Once done I'll let you know how it goes and then I can reward you the bounty. – alfonso Oct 12 '18 at 18:54
  • Your query involving ST_DWithin worked. The trouble I was having with the tabblock table was related to PostGIS. As of this moment, if you're using PostGIS 2.5, the tiger scripts are all messed up. Luckily, on my AWS instance I had PostGIS 2.4.4. The scripts are still messed up in 2.4.4, but you can easily fix it. For example, in KS, I had to edit the script by changing "ALTER TABLE tiger_staging.KS_tabblock" to "ALTER TABLE tiger_staging.KS_tabblock10". Seems like the census folks update their urls willy nilly and postgis breaks. Same applies to rest of the state and DC files. – alfonso Oct 22 '18 at 22:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.