[GIS] Finding all points in linestring that are contained by polygon using PostGIS

point-in-polygonpostgis

I have a linestring made up of 4D points and I need to determine which points in this linestring fall within a specific bounding box so that I can get the Z dimension of those points. I'm using Python 2.7 and Postgresql 9.3.4 with PostGIS 2.1.2. As of right now I'm first determining if the linstring overlaps the bounding box:

SELECT id FROM my_table WHERE ST_Intersects(my_geom, ST_GeometryFromText('POLYGON((minLon minLat, minLon maxLat, maxLon maxLat, maxLon minLat, minLon minLat))', 0))

Next, I'm dumping out all of the points of any linestring that intersects the bounding box:

SELECT ST_AsText((dp).geom) FROM (SELECT ST_DumpPoints(my_geom) as dp FROM my_table WHERE id like id_from_above) as foo

I then loop through the points and store any that are found inside the box (I'm pseudo-coding here):

for i in range(0, len(points)):
    SELECT ST_Contains(ST_GeometryFromText('POLYGON((minLon minLat, minLon maxLat, maxLon maxLat, maxLon minLat, minLon minLat))'), ST_GeometryFromText('points[i]'))

Last, I get the Z dimension of each of these stored points (again, pseudo-coding):

for i in range(0, len(stored_points)):
    SELECT ST_Z(ST_GeometryFromText('stored_points[i]'))

This method seems to work, however each linestring contains on the order of 25000 points. Since there are potentially many linestrings that will need to be checked, I'd like something that runs a bit faster.

Best Answer

Primary reason for the slower operation is due to unnecessary text conversions and out of sql processing. An optimised version of your own code is given below which should give you considerable speed enhancements. I tried to make it easier for you to understand the code. You should be able to run it directly from the sql console if your previous code was properly working.

SELECT ST_Z(pt) FROM 
(
  SELECT (dp).geom as pt FROM 
  (
    SELECT ST_DumpPoints(my_geom) as dp 
    FROM my_table 
    WHERE ST_Intersects(my_geom, ST_GeometryFromText('POLYGON((minLon minLat, minLon maxLat, maxLon maxLat, maxLon minLat, minLon minLat))', 0))
  ) as foo
  WHERE ST_Contains(ST_GeometryFromText('POLYGON((minLon minLat, minLon maxLat, maxLon maxLat, maxLon minLat, minLon minLat))'), (dp).geom)
) as foobar
Related Question