Sunday 30 October 2016

Postgis Self-intersection on geometryraster comparison


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?





No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...