20

I have a table with a column the_geom which contains data similar to:

0103000020E61000000100000005000000CE473AACFA071E40F27FB23340744740336FE841C6231E40873BED903F744740FC150A0ACE231E40D19E2684637647409C9B443D00081E409A9AF82664764740CE473AACFA071E40F27FB23340744740

Which when applying the function ST_AsEWKT(the_geom) returns:

SRID=4326;POLYGON((7.5077921782085 46.9082092877942,7.53493597966353 46.9081898840296,7.53496566473541 46.9249119938446,7.50781341296434 46.9249314035307,7.5077921782085 46.9082092877942))

I need to select all data that is within a 30km radius of a specific lat/long point, for example:

  • lat = 46.8167
  • lng = 6.9333

However whenever I tried to use ST_Distance(), I always received values less than 1, and using ST_DWithin() always returned true.

21

Please Check out the following query for PostgreSQL to get data within certain distance. I hope it will help.

SELECT *
FROM your_table
WHERE ST_Distance_Sphere(the_geom, ST_MakePoint(your_lon,your_lat)) <= radius_mi * 1609.34
  • 1
    Was able to get it working with: SELECT * FROM myTable WHERE GeometryType(ST_Centroid(the_geom)) = 'POINT' AND ST_Distance_Sphere( ST_Point(ST_X(ST_Centroid(the_geom)), ST_Y(ST_Centroid(the_geom))), (ST_MakePoint(6.9333, 46.8167))) <= 18 * 1609.34 – dan2k3k4 Nov 11 '13 at 16:57
  • That's great :) – Farhat Abbas Nov 11 '13 at 20:48
  • 2
    Just for anyone else wondering, the 1609.34 figure is meters per mile, which is the base unit postgres is using. So to do kilometers, obviously multiply by 1000. – 1mike12 Jun 27 '18 at 22:01
  • 1
    Pedantic note: it's 1609.344 exactly (by definition) – barrycarter Aug 17 '18 at 16:32
6

It sounds like you're storing your geometry in a geometry column, not a geography column.
That is fine, but the function ST_Distance will return measurements in projection units instead of always meters. In your case (4326), that will be degrees.
Just using a buffer with ST_Within won't work either, as the ST_Buffer will measured with degrees too.

You can either convert your data to use geography instead of geometry, or you can convert your point into some projection that uses meters, buffer, then convert back to 4326 to see what is within it:

SELECT
    *
FROM <your data>
WHERE ST_Within(the_geom, 
                ST_Transform(ST_Buffer(ST_Transform(ST_SetSRID(ST_MakePoint(6.9333, 46.8167), 4326), 3857), 30000), 4326)) = 1

That projects the point into 3857, which is a projection popular with web maps. Then it buffers it by 30,000 meters then reprojects it back to 4326 before passing it to ST_Within.

  • Except a pseudo-Mercator is unreliable for distance, so unless the data is close to the equator, the results will be off, especially with a distance of 30km. – Vince Jul 11 '17 at 1:49
6

In my world, using a custom SRID (for Google Maps) something like this worked:

SELECT * FROM addresses WHERE ST_DWithin(location, ST_SetSRID(ST_MakePoint(longitude, latitude), 3785), radius);

where the type of location is a geometry(Point,3785), and longitude, latitude, and radius are floats (e.g. -100, 44, 30 for 100W/44N/30 "units" -- see below)

See What is the best way to find all objects within a radius of another object? in the postgis docs:

The ST_DWithin(geometry, geometry, distance) function is a handy way of performing an indexed distance search. It works by creating a search rectangle large enough to enclose the distance radius, then performing an exact distance search on the indexed subset of results.

UPDATE: units are not miles for SRID 3785... they seem to be either radians or degrees or something like that. But the specification for my SRID says its units are either meters or degrees and it's definitely neither of those, at least not without some conversion:

alex=# select * from spatial_ref_sys where srid=3785; srid | auth_name | auth_srid | srtext | proj4text
3785 | EPSG | 3785 | PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]] | +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs

  • What's the difference between 3785 (in your post) and 3857? – tuned Jan 12 '16 at 17:35
  • They're different projections. 3875 vs 3857 -- I don't know if one is better than another – AlexChaffee Jan 13 '16 at 18:58
  • 1
    "EPSG 3785 has been deprecated in favor of the otherwise identical EPSG 3857" - github.com/rgeo/rgeo/pull/61 – Yarin Jan 25 '17 at 22:19
1

I think that this should work:

SELECT gid FROM table 
WHERE ST_DWithin(the_geom, ST_SetSRID(ST_Point(6.9333, 46.8167), 4326), 30000)
  • 2
    if you cast the_geom to geography that should work. st_dwithin(geography(the_geom),geography(<Point, 4326>), 30000) – cavila Sep 23 '15 at 21:11

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.