Querying Thousands of Points with ST_Value()?
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,
- Complex SQL, with simple library. "Complex" meaning container-types (rows, json, hstore, etc).
- 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
- Hstore, or
- 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,
We'll create points with this (it's just faster and more accurate),
ST_SetSRID(ST_MakePoint(lon, lat), 4326);
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,
- 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.
- 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