[raster] ST_Rescale changes extent

[strk]

According to documentation ( ​http://postgis.net/docs/RT_ST_Rescale.html ), ST_Rescale should only change the pixel scale, by retaining (or at most slightly enlarging) the extent of a raster.

While it mostly work that way, I've found this input to be handled differently, by retaining the same scale and halving the extent instead (!!)

0100000100000000000000F03F000000000000F0BF0000000000000000000000000000000000000000000000000000000000000000000000000A000A00640000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

Its metadata:

upperleftx | 0
upperlefty | 0
width      | 10
height     | 10
scalex     | 1
scaley     | -1
skewx      | 0
skewy      | 0
srid       | 0
numbands   | 1

Can be created with:

SELECT ST_AddBand(
 ST_MakeEmptyRaster(10, 10, 0, 0, 1, -1, 0, 0, 0),
 1, '8BUI', 0, 0
);                                                

[strk]

In other words, ST_Rescale does not change the scale, but just the size, although the parameters are floating point numbers. Noticeable are the "default geotransform" notices:

pg21=# with data as ( 
   SELECT ST_AddBand(
          ST_MakeEmptyRaster(10, 10, 0, 0, 1, -1, 0, 0, 0),
          1, '8BUI', 0, 0
        ) r 
 ), 
 meta as ( 
   SELECT st_metadata(st_rescale(r,2.0,-2.0)) m from data 
 ) select (m).* from meta;
NOTICE:  Raster has default geotransform. Adjusting metadata for use of GDAL Warp API
CONTEXT:  SQL function "st_rescale" statement 1
-[ RECORD 1 ]--
upperleftx | 0
upperlefty | 0
width      | 5
height     | 5
scalex     | 1
scaley     | -1
skewx      | 0
skewy      | 0
srid       | 0
numbands   | 1

[strk]

Another view on the matter, showing how origin 0,scaley*height is threated specially:

pg21=# select x, y, m.scalex, m.scaley from generate_series(-1, 1) x, generate_series(-1, 1) y, ST_AddBand(ST_MakeEmptyRaster(10, 10, x, y, 1, -1, 0, 0, 0), 1, '8BUI', 0, 0) r, st_metadata(st_rescale(r,2.0,-2.0)) m;
NOTICE:  Raster has default geotransform. Adjusting metadata for use of GDAL Warp API
CONTEXT:  SQL function "st_rescale" statement 1
 x  | y  | scalex | scaley 
----+----+--------+--------
 -1 | -1 |      2 |     -2
  0 | -1 |      2 |     -2
  1 | -1 |      2 |     -2
 -1 |  0 |      2 |     -2
  0 |  0 |      1 |     -1
  1 |  0 |      2 |     -2
 -1 |  1 |      2 |     -2
  0 |  1 |      2 |     -2
  1 |  1 |      2 |     -2
(9 rows)

[strk]

r11215 may be the cause of this

[strk]

I confirm the problem is with re-setting geotransform to the original one as the last step of warping when the default geotransform was changed to an arbitrary projection.

If a scale change was requested, the reset should take into account that change. This is likely to also affect rotation...

[strk]

I've a fix for this, but it breaks an existing test which I belive is broken, the one introduced here: http://trac.osgeo.org/postgis/changeset/11215

That test pretends that the output scale of a ST_Rescale with parameter 2,2 is 1, -1 as the input. Shouldn't it be 2, -2 instead ?

[strk]

I'm attaching the current patch I wrote to fix the reported case, but it needs more love

[strk]

r12928 fixes this in trunk. It does not make any attempt at handling translation or rotation of so called "default georeferenced rasters", is that supposed to be possible with _st_gdalwarp ? I'd guess so. It'd be better verified before backporting.

[strk]

Warping with a rotation still results in wrong output:

select x, y, m.width, m.height, m.skewx, m.skewy from generate_series(-1, 1) x, generate_series(-1, 1) y, ST_AddBand(ST_MakeEmptyRaster(10, 10, x, y, 1, -1, 0, 0, 0), 1, '8BUI', 0, 0) r, st_metadata(_st_gdalwarp(r,'NearestNeighbor', 0.125,NULL,NULL,NULL,NULL,NULL,1,1)) m;
 x  | y  | width | height | skewx | skewy 
----+----+-------+--------+-------+-------
 -1 | -1 |    10 |     15 |     1 |     1
  0 | -1 |    10 |     15 |     1 |     1
  1 | -1 |    10 |     15 |     1 |     1
 -1 |  0 |    10 |     15 |     1 |     1
  0 |  0 |     1 |     12 |     1 |     1
  1 |  0 |    10 |     15 |     1 |     1
 -1 |  1 |    10 |     15 |     1 |     1
  0 |  1 |    10 |     15 |     1 |     1
  1 |  1 |    10 |     15 |     1 |     1
(9 rows)

Trying to fix that I got the failures reported in #2912

[strk]

Ok #2912 has nothing to do with the changes in output skew computation, so I'm continuing with tests of the fix for that.

[strk]

For the record, this bug affects any raster having georeference of ul:0,0 scale:1,-1 skew:0,0 no matter having or not a SRID. I mention it because Bborie told on IRC that only SRID-less rasters were affected instead (that's not the case).

[strk]

NOTE: a webmercator tile for a raster of 1 meter per pixel may very well be common with upper-left corner 0,0

[strk]

I've noticed that ST_Width(rast) returns 10 and ST_Height(rast) returns 15, for all tiles (sounds good). So what is ST_Metadata.width and ST_Metadata.height ? Why are them different ?

[strk]

Please forget teh ST_Width vs. ST_Metadata.Width comment, it was my mistake, they do both return the same (bogus) answer for the warped raster.

[strk]

Well I backported r12928 to the 2.1 branch with r12951 and I'm closing this. Will file another one for the skew issue.