Monday 28 October 2019

PostGIS ST_Buffer Radius Help


I am trying to find all points within a five-mile radius of a given point. I have a query like this:


SELECT * FROM table WHERE ST_Contains(ST_Buffer(geomFromText('POINT(0 0)', 4326), ?), latlon)

I cannot figure out what I put in place of ? (radius) to get five miles. Everything is in EPSG 4326, and according to PostGIS documentation (as best as I can tell), my radius should be in meters. If I put in 12,070.0m (approximately 5-miles), I get matches half way across the country. Does anyone know what I am missing?



Answer



Because your data isn't projected - it is points on a spheroid - linear distances sort of don't make sense. Five miles at the equator is a much smaller angle than 5 miles on the arctic circle. But luckily PostGIS (>= 1.5) has the answer you're looking for:


SELECT * FROM table WHERE ST_DWithin(ST_GeogFromText('SRID=4326;POINT(0,0)'), geography(latlon), 12070);

It has a geography type that is designed for just this sort of thing. It's similar to geometry, but it only ever uses EPSG:4326, and there are far fewer functions that work with it.



In the sample above, I have called ST_GeogFromText() (There is also a ST_GeographyFromText(), and I'm not sure if there is a difference) on the point of interest (it might work with regular WKT, because the SRID parameter is redundant), and cast the latlon column to the geography type. If you're doing lots of these, it can be more efficient to create a geography column in your table and skip the cast entirely. Finally, ST_DWithin() can take geography parameters, and it does the right thing with linear distances.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...