[postgis-users] PostGIS Raster efficient enough to get height given coords from a large DEM?

[Bborie Park]

On Sun, Dec 4, 2011 at 6:55 AM, Stefan Keller <sfkeller at gmail.com> wrote:
> 2011/12/4 Bborie Park <bkpark at ucdavis.edu> wrote...
>>> Q2: Given the raster has a resolution (grid) of - say at best - 30m,
>>> does this function really interpolate, given the input is an arbitrary
>>> coordinate (I can't find a hint to that in the docs
>>> (http://postgis.refractions.net/documentation/manual-svn/RT_ST_Value.html))?
>>>
>>
>> There is no interpolation.  It literally is a pin stuck through the raster.
>
> You mean aka 'snap-to-grid', i.e. returning the closest grid point?
> Which grid point exactly?
>

Correct. It converts your geospatial value into a grid pixel x, y
(row, column) and then gets the value at the x, y.

> I now found a related thread here where elliott-22 asks "what if the
> yourLongitude, yourLatitude values provided in the query are not grid
> points in that they fall within a bounding pixel grid.  Only values
> for the intersection points are known so how does a non-intersecting
> point get assigned a value from the query?" (see
> http://old.nabble.com/Retrieving-Data-from-PostGIS-Raster-using-ST_Value-td32638248.html
>  at bottom of thread, Oct 17, 2011; 09:27pm)
>
> There Pierre answers:
>> Point outside the area covered by the grid are assigned NULL. There is not inter- or extra- polation done.
>
> But the question is about a query point (of float type) lying inside
> the perimeter but outside the grid nodes (probably integer based!
> These query points are inside the convex hull of the whole raster but
> they still don't "fall" on a regular grid point (see the figure of
> eliott-22 above).
>

You can't think that the value is only significant for the upper-left
corner of a pixel (regular grid point) as each pixel represents an
area (hence the geotransform values of scale, upperleft corner and
skew).  Since each pixel of a raster represents an area bounded by the
pixel's four corners with a singular value for that area (per band),
any point that falls within that pixel would have the same singular
value.

> What is exactly happening here in ST_Value()?
>

My response above.  ST_Value has the exact same behavior you'd expect
if you used the "Identify" tool in ArcGIS, ENVI, ERDAS, qgis and
clicked anywhere on a pixel.  It gets the value(s) (if more than one
band) of the pixel your cursor is over.

> Pierre also answered in the thread mentioned above my interpolation question:
>> Interpolating the elevation of a point is another story.
>> You could build a grid of 8, 15 or 24 points around your lat/long,
>> query the elevation at those points and use your algorithm to compute the interpolation.
>
> How do I really control which grid points are chosen - given an
> arbitrary (floating number) coordinate?
>

For what you need, you may want to convert an area of the raster
around your coordinates to geometries.

-- 
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
bkpark at ucdavis.edu