4

Currently I have an application using the psychopg2 Python library to query a database for elevation data. My current python implementation looks like the following:

   def GetElevation(lat, lon, cur):
    point = "'SRID=4326;POINT({} {})'::geometry".format(lat,lon)
    cur.execute("SELECT ST_Value(rast, {}) FROM dted0 WHERE ST_Intersects(rast, {});".format(point, point))
    return cur.fetchone()[0]

This works, but I am was curious if I could pass in an array of latitutdes and an array of longitudes. I tried building a "point string" with several thousand queries, but I get an error saying that I can only pass 100 arguments to the function.

Is there a way to pass all the points I need in a single query?

The reason I am doing this is because I assumed doing only one transaction with the database would be faster rather than looping for each query.

-------------------EDIT 1----------------------

In an attempt to implement the top rated post, I've done the following:

import json
json_arr = []
for lat, lon in zip(lats, lons): 
    json_arr.append({'lat:': lat, 'lon':lon})
json_str = json.dumps(json_arr)

Next I defined a function per his advice:

def GetElevationsBatch(json_elevations, cur):
    cur.execute( 
    """SELECT ST_Value(rast,z.point)
       FROM dted0
       JOIN (
          SELECT ST_SetSRID(ST_MakePoint(lat, lon), 4326) as point
          FROM json_to_recordset(%s) AS z(lon double precision, lat double precision)
       ) AS z
       ON ST_Intersects(rast,z.point)""", (json_elevations,))
    return cur.fetchall()

However, when I call the function I get I don't get good things:

a = GetElevationsBatch(json_str, cur)
print("a = ", a)
# Result: a =  None

-----------------------EDIT 2--------------------------

I've left the python plugin out for now until I can find the best query, so here is my latest attempt:

SELECT (ST_Dump(gv.geom)).geom, gv.val
FROM srtm , LATERAL ST_Intersection(rast, 'SRID=4326;MULTIPOINT(-111.305048568 38.0601633931,-111.822991286 38.6025320796,-111.136977796 38.3631992596,-111.206470006 38.971228396)'::geometry)  AS gv
WHERE ST_Intersects(rast, 'SRID=4326;MULTIPOINT(-111.305048568 38.0601633931,-111.822991286 38.6025320796,-111.136977796 38.3631992596,-111.206470006 38.971228396)'::geometry);

In the above query, I am attempting to use the notion of MULTIPOINT geometries. However, I found that these calls "work", but they are much slower than my original method of simply querying for every point. I understand that I am not running on any amazing hardware, but should calls really take on the order of seconds for a simple elevation query? Seems to me there is something awry here. It is taking nearly 3 seconds to retrieve 4 points with the above code. To compare to my original method, it only takes about 100ms to retrieve 4 points.

Shouldn't a solution in which I only query the database once be quick than one where I have to query several times?

  • Are you hoping to return a dataset of only elevations, or of elevations with their corresponding lat/long coordinates? – Nate Wanner Nov 21 '16 at 22:13
  • I only need elevations because I already know what the lat/longs are since I'm querying for them. – user985030 Nov 21 '16 at 22:50
  • You could ST_Dump a MultiPoint in your SQL, so they you're passing a MultiPoint but getting back lots of rows of z's. – Paul Ramsey Nov 21 '16 at 23:12
  • Given that you are trying to pass more than 100 points, is it really a set of points you are interested in or a specific area? Maybe if you state what you plan to do with these values once you get them, the job is easier. – LR1234567 Nov 22 '16 at 5:09
  • ST_MakePoint takes (lon, lat), not (lat, lon) – BradHards Dec 3 '16 at 23:10
1

Okay I thought of another reason why my original answer might be slow. If the bounding box of your point covers enough area, it would produce a lot of rasters that require checking the Slow way.

So here is another answer, still using multi-point but doing your original single check approach:

SELECT dp.geom, ST_Value(dted0.rast, geom) AS val
   FROM ST_Dump(your_multi_point_here) AS dp 
     JOIN  dted0 ON ST_Intersects(dted0.rast,dp.geom) ;
|improve this answer|||||
  • This is much faster than your previous solution. With this solution I was able to do about 10,000 elevations per second running postGIS in a Docker container. – user985030 Dec 5 '16 at 23:03
2

Yes. You can. First. Don't use .format() and Python curly brace syntax. Use the Psycopg placeholders. In the docs.

Warning Never, never, NEVER use Python string concatenation (+) or string parameters interpolation (%) to pass variables to a SQL query string. Not even at gunpoint.

Second, you need to solve the problem of how to represent multiple lat/long. You can do this various ways, two popular methods,

  1. Complex SQL, with simple library. "Complex" meaning container-types (rows, json, hstore, etc).
  2. Simple SQL, with complex library.

Pyscog2 is a simple library. It provides very little abstraction over sql. Perl's SQL::Abstract is a bit more complex, and ORM's are comparatively deep voodoo magic. So with only a simple library like pyscog, your options are to have it serialize the options into

  1. Hstore, or
  2. JSON.

Let's look at hstore. Hstore is the default serialization type for Python's dict. It's not ideal: if you have a dict of lat=>lon, what will you do if two lats are the same? So, we have to use JSON, which supports an array.

Next lets draw up a method,

  1. We'll create points with this (it's just faster and more accurate),

    ST_SetSRID(ST_MakePoint(lon, lat), 4326);

  2. We'll get the data in with json_to_recordset. With this we just need to send in a json array, '[{ "lat": float, "long": float }...]'

Now we just need to do something like this....

SELECT ST_Value(rast,z.point)
FROM dted0
JOIN (
  SELECT ST_SetSRID(ST_MakePoint(lat, long), 4326) as point
  FROM json_to_recordset(%s) AS z(long double precision, lat double precision)
) AS z
  ON ST_Intersects(rast,z.point)

Some of the people are pointing out you can do ST_Intersection instead that is true. Let's review,

  1. ST_Intersection — (T) Returns a geometry that represents the shared portion of geomA and geomB. The geography implementation does a transform to geometry to do the intersection and then transform back to WGS84.
  2. ST_Intersects — Returns TRUE if the Geometries/Geography "spatially intersect in 2D" - (share any portion of space) and FALSE if they don't (they are Disjoint). For geography -- tolerance is 0.00001 meters (so any points that close are considered to intersect)

So we've already solved the major problem of getting the lat long coordinates into the database. We solved this problem by serializing those coordinators into JSON. It's important to note that this problem could have also been solved by passing in,

ST_GeomFromEWKT($$SRID=4326;MULTIPOINT (10 40, 40 30, 20 20, 30 10)$$);

To use this method with ST_DumpValues the query would look something like

SELECT *
FROM dted0
FULL OUTER JOIN (
  SELECT rast,*
  FROM ST_DumpValues(
    ST_Intersection(rast, ST_GeomFromEKWT(%s))
    , band
  )
) AS z
  ON z.rast = rast;

This should return something like rast|band|valarray

|improve this answer|||||
  • 1
    Have a typo, I think you meant the ON to be ST_Intersects and not ST_Intersection. – LR1234567 Nov 22 '16 at 6:42
  • I think your ST_Intersection should be ST_Clip. – LR1234567 Nov 23 '16 at 5:42
  • I've tried both of your solutions and neither one of them works. Most recently, I tried your second solution and I get an error "rast does not exist". Your second reference to it in ST_DumpValues does not work. – user985030 Dec 3 '16 at 19:44
1

I'm thinking your best bet is using ST_Intersection. As Paul mentioned, you can wrap all your points in a multipoint if you are really looking at individual points. If its an area you are really looking for the dem for, this will work for that too.

http://postgis.net/docs/manual-2.3/RT_ST_Intersection.html

SELECT gv.geom, gv.val
   FROM dted0 , LATERAL ST_Intersection(rast, {})  AS gv
   WHERE ST_Intersects(rast, {});

Your {} can be a MULTIPOINT, LINESTRING, POLYGON, MULTI.. any valid 2D geometry.

Note that if your geometry has lots of different values, the geometry will be split across pixel values. So if you pass in a MULTIPOINT, you might get back a set of POINTS and MULTIPOINTS.

If you really want them to be Points, then employ Paul's idea of ST_Dump Like so:

SELECT (ST_Dump(gv.geom)).geom, gv.val
   FROM dted0 , LATERAL ST_Intersection(rast, {})  AS gv
   WHERE ST_Intersects(rast, {});
|improve this answer|||||
  • ST_Intersection(raster, geometry) returns a raster. ST_Dump takes a geometry. I think you need ST_DumpValues() – Evan Carroll Nov 22 '16 at 15:25
  • 1
    No ST_Intersection(raster,geometry) returns a set of geomvals. I think you are thinking about ST_Intersection(raster,raster). A geomval is a composite object consisting of a geometry and a pixel value. – LR1234567 Nov 22 '16 at 18:27
  • I found this solution to be very slow. When I ran a quick comparison of methods, my original function could query a single point in about 20ms, your multipoint solution takes about 1s. I then attempted a single call with 5 lat,lon pairs and this took about 3s. It seems that this solution is much slower than just querying everyone point individually. – user985030 Dec 3 '16 at 22:25
  • Can you provide the query plan of this. Puzzled why that is much slower. Might be index is not being used for some reason. Can you also try flipping the order of args to ST_Intersection(geom, rast). As I recall I think one does the operation in raster space and one in geometry or used to. So that might be causing some difference depending on size of rasters. – LR1234567 Dec 5 '16 at 17:02
  • Not sure by what you mean, "query plan"? – user985030 Dec 5 '16 at 23:04
0

I think there will be another way to query multiple points on raster datasets.

WITH pairs(x,y) AS (
VALUES
(70.5084, 45.1985),
(70.5094, 45.1995),
(70.5084, 45.1975),
(70.5074, 45.1965),
(70.5064, 45.1955),
(70.5054, 45.1945)
) 
 SELECT 
    ST_Value(rast, ST_SetSRID(ST_MakePoint(x, y), 4326)) AS height
 FROM dem_raster rs
    CROSS JOIN pairs
    WHERE ST_Intersects(rs.rast,  ST_SetSRID(ST_MakePoint(x, y), 4326));

My DEM data will respond this query within < 50 ms.

|improve this answer|||||

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.