Postgis SQL for nearest neighbors

Question

I'm trying to calculate the nearest neighbors. For that, I need to pass a parameter to limit the maximum distance from the neighbors. For example, which are the nearest neighbors within a radius of 1000 meters?

I did the following :

I created my table with the data:

id | name | latitude | longitude

After that, I executed the following query :

SELECT AddGeometryColumn ( 'public' , ' green ', ' geom ' , 4326 , ' POINT' , 2 );

UPDATE season
SET geom = ST_Transform(ST_PointFromText ('POINT (' || longitude || ' ' || latitude || ')', 4269), 4326);

First question, Is the SRID of Brazil 4326? What would be 4269 ?

Second question, by doing the following SQL

SELECT id, name
FROM season
WHERE ST_DWithin (
                 geom ,
                 ST_GeomFromText ('POINT(-49.2653819 -25.4244287 )', 4326),
                 1000
                 );

This returns nothing. From what I understand, this SQL would further point the radius of the maximum distance, right?

It appears if you put 1000 results for 100000000, all my entries appear .

So, I wonder what is wrong here?

Answer 1

First, If you are using latitude, longitude, you need to use 4326.

UPDATE season SET geom = ST_PointFromText ('POINT(' || longitude || ' ' || latitude || ')' , 4326 ) ;

Then you create an index on the geom field

CREATE INDEX [indexname] ON [tablename] USING GIST ( [geometryfield] ); 

Then you get the kNN neightbors:

SELECT *,ST_Distance(geom,'SRID=4326;POINT(newLon newLat)'::geometry) 
FROM yourDbTable
ORDER BY
yourDbTable.geom <->'SRID=4326;POINT(newLon newLat)'::geometry
LIMIT 10;

This query will take advantage of kNN functionality of the gist index (http://workshops.boundlessgeo.com/postgis-intro/knn.html).

Still the distance returned will be in degrees, not meters (projection 4326 uses degrees).

To fix this:

SELECT *,ST_Distance(geography(geom),ST_GeographyFromText('POINT(newLon newLat)') 
FROM yourDbTable
ORDER BY
yourDbTable.geom <->'SRID=4326;POINT(newLon newLat)'::geometry
LIMIT 10;

When you calculate the ST_distance use the geography type. There distance is always in meters:

http://workshops.boundlessgeo.com/postgis-intro/geography.html

All this functionality will probably need a recent Postgis version (2.0+). I am not sure though.

Check this for reference https://gis.stackexchange.com/questions/91765/improve-speed-of-postgis-nearest-neighbor-query/

Answer 2

If you want to find information in Postgis on spatial reference systems look in the table spatial_ref_sys. Based on that and to answer you first question,

select * from spatial_ref_sys where srid=4629

returns the following:

4269 | EPSG      |      4269 | GEOGCS["NAD83",DATUM["North_American_Datum_1983"
,SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["
EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01
745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]] | +proj=long
lat +ellps=GRS80 +datum=NAD83 +no_defs

This indicates that it is a North American Datum, and not Brasil, using the GRS80 ellipse, which is not the same as the WGS84 one used by 4326, which is the standard used by GPS. This information is what is used by the Proj4 projections library, which is included as part of a Postgis build on Postgres.

I searched for Brasil in spatial_ref_sys,

select * from spatial_ref_sys where auth_name ilike '%Bras%';

and this returned no results.

I did find a reference to 29101 on the web, which uses a South American datum, with a Brazilian Polyconic projection and units in metres.

To convert you original example to this, you could use:

select ST_Astext(ST_Transform(ST_SetSrid(ST_Makepoint(-49.265, -25.424),4326), 29101));

which returns:

POINT(5476247.04359163 7178517.77380949)

and would allow you to use ST_Distance or ST_DWithin in meters. I only put the ST_AsText so you can see the output, you obviously wouldn't actually use that in an update using ST_Transform.

Note the use of ST_Makepoint instead of concatenating coordinates -- the same thing, but easier to read, imho.

To answer your second question, because you are effectively working in 4326, lat/lon, which ranges from -180 to 180 and -90 to 90, if you used 1000 as the distance in ST_DWithin, you will get every record returned, because the 1000 is not meters, but degrees, so you would have to use something like the system above.

If you want to take distance in lat/lon, have a look at the haversine formula which will work fine for short distances.