Monday 22 February 2016

postgis - How to use ST_Intersection?


Here's a quick summary about what I'm trying to do: I have 3 tables in Postgres, 'a' and 'b', each have a Polygon column, and 'c' has a Point column. What I'm trying to do here is to get the geometries intersections between 'a', 'b' and 'c', and to display such geometries on an OpenLayers vector layer.


I already know how to display any kind of geometry from a String in OpenLayers, but I'm having troubles with the PostGIS' ST_Intersection function, I'm doing this:


SELECT ST_Intersection(a.geom, b.geom) as inter from a, b;

where a.geom and b.geom are both the geometry columns, and I get this error message:


NOTICE:  TopologyException: found non-noded intersection between 515172 2.14408e+06, 497067 2.13373e+06 and 501321 2.13546e+06, 471202 2.14843e+06 500621 2.13576e+06 

ERROR: GEOS Intersection() threw an error!

Also I tried to express the resultant geometry as text using ST_AsText like this:


SELECT ST_AsText(ST_Intersection(a.geom, b.geom)) as inter from a, b;

but it send me this error message:


HINT: No function matches the given name and argument types. You might need to add explicit type casts.

I don't know what I'm doing wrong, I just want to get the Polygons' WKT to display it on OpenLayers, here's how I display a geometry from a WKT:


                    var in_options = {

'internalProjection': new OpenLayers.Projection("EPSG:4326"),
'externalProjection': new OpenLayers.Projection("EPSG:4326")
};

var fea= new OpenLayers.Format.WKT(in_options).read(data); //data is the string with the WKT
vectorLayer.addFeatures([fea]); //this piece of code works great
map.zoomToExtent(bounds);

UPDATE: I tried the next:


SELECT ST_Intersection(a.geom, b.geom) as intersect_ab FROM a INNER JOIN b ON 

ST_Intersection(a,b) WHERE ST_Overlaps(a.geom, b.geom)
AND ST_isvalid(a.geom)='t' AND ST_isvalid(b.geom)='t';

but I get the next error message:


ERROR: Function st_intersection(a,b) does not exist.
HINT: No function matches the given name and argument types. You might need to add explicit type casts.

I added the isvalid to verify only valid polygons are being evaluated, but it's telling the error is in the ST_Intersection(a,b), both a, b and c have the same SRID so I'm really confused, sorry if I'm asking too much, but I'm quite new with PostGIS so I hope I'm not bothering you a lot. Thanks.



Answer



My guess would be that it fails if the intersection returns NULL. So you should add a where clause checking if there actually is an intersection before you try to create the WKT.



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