Monday, 24 September 2018

postgresql - Get geometry from coordinates


I have a table with geometry column called point, saved in SRID 4326. I want to get the geometry by its coordinates, but something didn't work. I suspect database configuration or something.


Here is what I do


select point from geom_data limit 1

And I get "0101000060E610000000000020558A5EC0000000C0F3B8424000B0CFCFA4717342"


This as text is



select AsText(point) from geom_data where point = '0101000060E610000000000020558A5EC0000000C0F3B8424000B0CFCFA4717342'

And the result is "POINT(-122.161445617676 37.444938659668)"


Trying this


select * from geom_data where ST_Contains(point, ST_GeomFromEWKT('SRID=4326;POINT(-122.161445617676 37.444938659668)'))

Didn't give a result. Tried this


select ST_PointFromText('POINT(-122.161445617676 37.444938659668)', 4326)

and get "0101000020E61000000F000020558A5EC0040000C0F3B84240" which is different.



Selecting within bounding box didn't worked either


select * from geom_data where
st_x(point) >= '-122.161445617670' AND st_x(point) <= '-122.161445617679'
AND st_y(point) >= '37.444938659660' AND st_y(point) <= '37.444938659669'

Any ideas?



Answer



There are two parts to your question really. First, the binary discrepancy between the point in the database and the point you supply, and second, trying to do a spatial query on the points.


To address the discrepancy problem, it turns out that the point stored in your database is a 2D+M point. Running this query:


SELECT ST_AsEWKT('0101000060E610000000000020558A5EC0000000C0F3B8424000B0CFCFA4717342');


Yields:


SRID=4326;POINTM(-122.161445617676 37.444938659668 1336176082171)

This is fine for spatial queries, because PostGIS will ignore the M parameter, which is used as a general per-vertex floating point attribute. Originally it stood for "measure", but it can represent anything you fancy.


The second problem can be solved in two ways. You can use the ST_Equals() function:


SELECT ST_Equals(ST_GeomFromEWKT('SRID=4326;POINTM(-122.161445617676 37.444938659668 1336176082171)'), ST_PointFromText('POINT(-122.161445617676 37.444938659668)', 4326));

or


SELECT ST_Equals(ST_AsEWKT('0101000060E610000000000020558A5EC0000000C0F3B8424000B0CFCFA4717342'), ST_PointFromText('POINT(-122.161445617676 37.444938659668)', 4326));


Both of which return true.


Or use ST_Within() if there is some tolerance needed:


SELECT ST_DWithin(ST_GeogFromText('SRID=4326;POINTM(-122.161445617676 37.444938659668 1336176082171)'), ST_GeogFromText('SRID=4326;POINT(-122.161445617676 37.444938659668)'), 1.0)

Notice here I'm using the geography type, this is so you can specify the distance in metres on unprojected data. If you just used the geometry type, you'd have to specify the tolerance in degrees, which is fraught and generally of not much use, or project the data first which is potentially inefficient.


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