3

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?

1
  • I did not find a solution to getting rid of the notices, so I decided to ignore them. Using C++ this can be done by setting an "empty function" as dummy Notice Processor using "PQsetNoticeProcessor(conn, dummy, 0);".
    – Niklas
    Jan 27, 2017 at 17:47

1 Answer 1

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
2
  • Thank you for your answer. I only just now found time to continue this project. How can I find out if my data is already tiled? Is it possible to tile existing raster data or do I need to reimport it? Just adding an Index did not seem to improve speed.
    – Niklas
    Jan 27, 2017 at 17:41
  • 1
    For the record: Indexing gives a huge performance boost. After a quick 'n dirty evaluation using ASTER GDEM on Baden-Württemberg (county in southern Germany) tiling with 10x10 and indexing seems to be working best. Even smaller tiles seem to have similar performance but take far more disk space.
    – Niklas
    Feb 14, 2017 at 9:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.