I'm a Rails developer with minimal GIS experience, working with PostGIS and the RGeo gem.
I'd like to find all of the cows within 500 meters of my current location. So, I'll create a model for my cows:
class Cow < ActiveRecord::Base
set_rgeo_factory_for_column(:latlon, RGeo::Geographic.spherical_factory(:srid => 4326))
I'll then seed my database with enough cows to keep me busy:
5000.times{ Cow.new(latlon:"POINT(#{40.735767 + rand() * 0.1 } #{-77.919674 + rand() * 0.1 })" ).save }
To find the cows within 500 meters of my location, I'll use the ST_DWithin PostGIS function that all the tutorials suggest:
cows_found = Cow.find_by_sql("SELECT * FROM cows WHERE ST_DWithin(latlon, '#{search_latlon}', 500)")
The unhappy result is the following:
Moving to the East or the West, the 500 meter radius threshold appears to be observed correctly, but moving to the North or the South, I am finding cows as far as two kilometers away.
Now, obviously I'm extremely naive in my approach, but I'm not completely clueless; I understand I'm populating my cows in a spheroidal projection and handing the coordinates returned from ST_DWithin to Google's Mercator. I expect distortion. However, this is larger than I anticipated and suggests to me that I've fundamentally misunderstood something.
I had thought that in my use of ST_DWithin I was passing as arguments two geographies, and that distance would therefore be specified in meters. I was able to set the use_spheroid boolean argument to false without PostGIS complaining, which I do not think I should be able to do if I were in fact comparing geometries. In fact, if I use psql to look at the latlon column, it is of type geography:
latlon | geography(Point,4326) |
However, perusing related questions, I see several answers that state that using ST_DWithin with SRID 4326 results in the distance threshold being specified in radians. (Edit: I now understand this would only be the case if one was using 4326 geometries, not 4326 geographies.)
So, given the goal of locating cows within a given number of meters, and tolerating an error of tens of meters (nobody will die if we find a few extra cows) what would be the best way to proceed?
Throw out the 4326 spatial reference system, and use a different SRID?
Use 4326, but calculate a conversion from meters to radians, for use with ST_DWithin? (This conversion would obviously be subject to the specific latitude and longitude from which I am searching.)
Any suggestions would be appreciated.
Edit:
Both arguments to the ST_DWithin call are 4326 geographies. The first argument, latlon, is a 4326 geography column, and the second argument is drawn from another column of the same type. A commenter's (very perceptive) suggestion that I ensure the second argument was in fact a geography by using ST_GeogFromText led to precisely the same end result.
So, I have high degree of confidence that I am using the correct version of ST_DWithin with the correct parameters. As a sanity test, I repeated exactly the same procedure (generating a dense set of latlons and performing ST_DWithin upon it) very near the equator. This is the result:
As you can see, the points returned by ST_DWithin are much less distorted. More importantly, they are all within 500 meters of the central point. The results that I am able to obtain at the equator are precisely in line with my expectations. However, at roughly the 41st parallel North, I am getting the distortion observed in the first example.
Now, I expect distortion from handing my 4326 coords to Google Maps (who internally transform them to their own web Mercator) and I expect some distortion from using 4326 (PostGIS In Action warns that local distance measurements will be 'bad'). I expect these distortions to become more severe as I move away from the equator. The first map example in this question is obviously subject to both of those distortions.
So my question is: Is the level of distortion I am seeing in the first map example actually a correct result? (Meaning, the result one would expect given the limitations of the 4326 spheroid model and distortion from the projection to web Mercator?)
If this is a correct result, I would be eager to hear suggestions regarding a PostGIS strategy to retrieve points within a given radius that was less subject to these distortions, or a strategy to ameliorate this distortion.
