1

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.

3

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
0

Something like this.

WITH poly AS (
  SELECT ST_MakeEnvelope(...) AS geom
),
pts AS (
  SELECT (ST_DumpPoints(l.geom)).geom AS geom, id AS line_id
  FROM poly p JOIN lines l ON ST_Intersects(p.geom, l.geom)
)
SELECT line_id, ST_Z(pts.geom)
FROM pts JOIN poly ON ST_Contains(poly.geom, pts.geom)
  • 1
    Thanks for your response. Can you elaborate a bit for me? As in what your variables are in relation to mine. Specifically I'm having trouble with the JOIN lines l part. – abe678 Jun 4 '14 at 13:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.