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

[Stefan Keller]

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?

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

What is exactly happening here in ST_Value()?

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?

Stefan


2011/12/4 Bborie Park <bkpark at ucdavis.edu>:
> On Sat, Dec 3, 2011 at 7:16 PM, Stefan Keller <sfkeller at gmail.com> wrote:
>> Hi Bborie
>>
>> Many thanks for your response.
>> To summarize: My goal is to get the height (in meters) given an
>> arbitrary lat/lon coordinate from the latest free digital elevation
>> model (like mentioned here
>> http://slashgeo.org/2011/12/01/New-Free-Elevation-Dataset-named-Global-Multi-resolution-Terrain-Elevation-Data-2010-GMTE
>> ).
>>
>> The approach using PostGIS to access raster data is
>> 1st. to convert an ASCII GRID text file (an ESRI format) to GeoTIFF
>> (not tiled nor "pyramide"-ized).
>> Then (2nd) to load it into PostGIS, like described here
>> http://trac.osgeo.org/postgis/wiki/WKTRasterTutorial01 .
>>
>
> I don't remember but I think GDAL supports ASCII GRID so you don't
> need an intermediary file.  The loader uses GDAL to convert them into
> the PostGIS Raster format.
>
>> I still have some basic questions about this interesting approach:
>>
>> Regarding my question of "rules of thumb" in order to set the tile
>> size you answered:
>>> Generally speaking, the smaller the tile size, the faster the operation.  In
>>> my testing of my rasters, 200x200 is slower than 100x100 which is slower
>>> than 50x50.  The smaller tile size results in more rows in the table BUT the
>>> performance dramatically increases.  I can't remember of the top of my head
>>> but the performance gain was linear or better.
>>
>> Thats a good hint.
>> Q1: But are there no calculcation rules on how to estimate an ideal
>> absolute size (in pixels)?
>>
>
> Nope.  Testing is the only way as there are other possible factors
> involved.  And even my estimates above are only guidelines.  I know
> I'm going to have to test for my production database servers before I
> know what works best for them.  In a pinch, I'd say you'd be happy
> between 30 - 70.
>
>> Regarding a query given lat/lon you proposed:
>>> Assuming you test to see what the ideal (or "good-enough") tile size
>>> is for the elevation rasterset, you should see good performance using
>>>
>>> ST_Value(rast, ST_SetSRID(ST_MakePoint(long, lat), srid))
>>>
>>> You could eliminate the ST_SetSRID and ST_MakePoint by passing a
>>> correct EWKT instead
>>>
>>> ST_Value(rast, 'YOUREWKT'::geometry)
>>
>> 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.
>
>> Q3: If yes, I assume its linear. How could I proceed to get a result
>> using non-linear interpolation (like Kriging)?
>>
>
> I don't know if you can using only PostGIS raster.  If you're able to
> use PL/R as well, you should be able to tap into their kriging tools.
>
> -bborie
>
> --
> Bborie Park
> Programmer
> Center for Vectorborne Diseases
> UC Davis
> 530-752-8380
> bkpark at ucdavis.edu
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users