Querying Thousands of Points with ST_Value()?

Question

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?

Answer 1

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

Answer 2

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) ;

Answer 3

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, {});

Answer 4

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.