[raster] Coverage topology broken by ST_Rescale

[strk]

As described here: ​http://lists.osgeo.org/pipermail/postgis-devel/2014-July/024337.html

The ST_Rescale function breaks coverage topology introducing gaps between a tile and another whereas ST_Resize does not.

The test dataset is here: ​http://www.naturalearthdata.com/downloads/10m-raster-data/10m-cross-blend-hypso/

Imported with this command:

raster2pgsql -C -I -s 4326 -Y -t '256x256' HYP_LR.tif hyp_lr

It looks like a bug for ST_Rescale to not maintain topology as the manual says that the tile extent could be slightly modified to accomodate for subpixel scaling.

I hadn't tested with older versions

[strk]

See also #2800

[strk]

Oops, I meant see #2833 (2800 is this one)

[strk]

May be related to #2251

[strk]

Adding some debugging lines I find that ST_Resize (the working one, being passed percent width/height) ends up calling rt_raster_gdal_warp with an integer width/height like this:

NOTICE:  rt_raster_gdal_warp user-defined width 3, computed xscale 0.6144
NOTICE:  rt_raster_gdal_warp user-defined height 3, computed yscale 0.6144

While ST_Rescale ends up calling rt_raster_gdal_warp with a float scale like this:

NOTICE:  rt_raster_gdal_warp user-defined scale x=0.6, y=0.6
NOTICE:  rt_raster_gdal_warp computed dims w=3, h=3

The source scale of the raster is: 0.00719999999999928 (both X and Y). The source size is 256 x 256.

ST_Resize ends up being called with percent parameters 0.0119999999999988, being the result of source scale / target scale (0.6). The full query is:

SELECT ST_Resize("rast", CASE WHEN abs(ST_ScaleX("rast"))::float8 < 1 THEN abs(ST_ScaleX("rast"))::float8/1 ELSE 1.0 END, CASE WHEN abs(ST_ScaleY("rast"))::float8 < 1 THEN abs(ST_ScaleY("rast"))::float8/1 ELSE 1.0 END) AS geom,"rid" FROM dataraster

ST_Rescale gets called with 0.6 instead.

The result from ST_Resize retains the original extent while the result from ST_Rescale makes it slightly smaller:

with r as ( select rast from hyp_1250m limit 1 ) select box2d(st_envelope(rast)) as orig, box2d(st_envelope(st_resize(rast,.0119999999999988,.0119999999999988))) as resized, box2d(st_envelope(st_rescale(rast,0.6,0.6))) as rescaled from r;
-[ RECORD 1 ]-------------------------------------
orig     | BOX(-180 88.1568000000002,-178.1568 90)
resized  | BOX(-180 88.1568000000002,-178.1568 90)
rescaled | BOX(-180 88.2,-178.2 90)

I guess all of this is expected, which would make the confusion simply due to lack of documentation. It should be stated in the manual page when ST_Resize retains extent (I guess when passed a factor of original size, but not when passed fixed size) and when (if ever) ST_Rescale retains extent.

[strk]

What strikes me is that ST_Rescale takes less time than ST_Resize.

explain analyze select st_rescale(rast,.6,.6) from hyp_1250m; -- rows=19208
 Total runtime: 29388.881 ms
explain analyze select st_resize(rast,.0119999999999988,.0119999999999988) from hyp_1250m; -- rows=19208
 Total runtime: 37063.069 ms

[strk]

Now I see that the additional time taken by ST_Resize is due to the function being implemented in in pl/pgsql and invoking ST_Metadata to fetch original width/height of the raster, and computing the new one.

Indeed the version taking integer width/height takes the same time of st_rescale. Closing this one as invalid, and let's refer to #2833 for improving the documentation

[strk]

I'm actually reopening this because the manual page for ST_Rescale says:

When the new scalex or scaley is not a divisor of the raster width or height, the extent of the resulting raster is expanded to encompass the extent of the provided raster

So the extent of the raster coming out from ST_Rescale should be _larger_ rather than _smaller_ than the input raster, as it's happening now.

[robe]

Milestone renamed