Postgis Self-intersection on geometry<->raster comparison

Question

I am rather new to postgis and try to get a (weighted) average height of an area around a point. I am running Postgresql 9.5.5 with Postgis 2.2.1 on Ubuntu. I have imported ASTER GDEM elevation information as a raster and built the following functions, following this tutorial:

WITH statics AS (
    SELECT (ST_DISTANCE(CAST (ST_Project(CAST(ST_SETSRID($1, 4326) AS geography), $2, radians(90.0)) AS geometry), CAST(ST_SETSRID($1, 4326) AS geometry))) AS exp_fact
    -- calculate degrees to expand the point with, $2 is given in meters
)
SELECT SUM(sub.height * sub.weight) / SUM(sub.weight)
FROM (
    SELECT height, ((1 - (dist / $2))) as weight
    FROM get_elev_points_within_geom(ST_MAKEVALID(ST_BUFFER(ST_SETSRID($1, 4326), (SELECT exp_fact from statics))), $1)
    WHERE dist <= $2
) as sub

where get_elev_points_within_geom(geo geometry, center geometry) is:

SELECT ST_Y(ST_CENTROID((a.neighborhood).geom)),
    ST_X(ST_CENTROID((a.neighborhood).geom)),
    (a.neighborhood).val,
    ST_DISTANCE(CAST (ST_CENTROID((a.neighborhood).geom) AS geography), CAST ($2 AS geography)) as dist_m
FROM(
    SELECT ST_INTERSECTION(elevation.rast, $1) AS neighborhood
    FROM elevation
    WHERE ST_INTERSECTS(elevation.rast, $1)
) AS a
ORDER BY dist_m
LIMIT 1000;

Running e.g. this query:

SELECT * FROM get_elev_average_weighted_by_dist_at(ST_MAKEPOINT(8.4305360482003504, 49.013877361340597), 50);

I get:

NOTICE:  Ring Self-intersection at or near point 8.4315277777777791 49.016527777777782
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
SQL function "get_elev_points_within_geom" statement 1
SQL function "get_elev_average_weighted_by_dist_at" statement 1
Total query runtime: 145 msec
1 row retrieved.

I have traced the Problem down to to:

SELECT (ST_INTERSECTION(elevation.rast, ST_SETSRID(ST_GeomFromText('POLYGON((8.4298 49.0131,8.4298 49.0145,8.4312 49.0145,8.4312 49.0131,8.4298 49.0131))'),4326))).geom
FROM elevation
WHERE ST_INTERSECTS(elevation.rast, ST_SETSRID(ST_GeomFromText('POLYGON((8.4298 49.0131,8.4298 49.0145,8.4312 49.0145,8.4312 49.0131,8.4298 49.0131))'),4326))
LIMIT 14;

which will return with the same error message (but different trace, obviously). Interestingly, if I reduce the limit to 13, I will not get the Message.

Also this will run for about 150ms per Query on my machine. This is quite annoying if I query for lots of points.

How can I improve on efficiency?

Answer 1

In regards to performance, maybe this is helpful for you: http://www.geonet.ch/postgres-postgis-of-rasters-and-geojsons/ The quintessence of this blog post in regards to performance is, that the best performance is achieved by tiling and indexing raster databases in PostGIS:

Test Raster Extent: 5000 X 2000 = 10’000’000 cells

• Querying the raster no tiling, no index: 522ms
• Querying the raster with tiling, no index: 227ms
• Querying the raster with tiling, with index: 84ms