Sunday, 22 February 2015

postgis - How to add X/Y coordinates and count of overlapping points in a postgresql table?


How to add to my postgres/postgis table the following columns



X coordinate | Y coordinate | count (number of points in same location)


from an existing table populated with geometries type POINT



Answer



The "same location" part of your question is problematic, since locations are stored as floating point values, and due to computer operating system rounding errors you will get two points that you consider at the same location with very slightly different X,Y.


The correct way to do this is to define a buffer distance within which you consider two points as equal. The PostGIS function that does this is ST_Dwithin(). So to complete Hicham Zouarhi's comment/answer (assuming your table is called "pts" and it has a column "the_geom", and a primary key column "id":


ALTER TABLE pts ADD COLUMN x DOUBLE PRECISION, y DOUBLE PRECISION, cnt INTEGER;
UPDATE pts SET x=ST_X(the_geom), y=ST_Y(the_geom);
UPDATE pts SET cnt=(SELECT count(*) FROM pts AS p1 JOIN pts AS p2 ON ST_Dwithin(p1.the_geom, p2.the_geom,0.001) WHERE p1.id>p2.id);

I chose a distance for ST_Dwithin of 0.001. If your data is in a meter based projection, then the buffer will be 1 mm. However, If your data is in a geographic, long/lat based CRS, then this buffer will be 0.001 degrees, or about 100 meters. So in a degree based CRS, you probably will want a much smaller distance.



The WHERE part of the UPDATE is to insure that you don't get double counts: i.e. point 1 is within 1mm of point 2, and point 2 is within 1 mm of point 1. By comparing the id values, you avoid this duplicate count.


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...