[raster] Problem with ST_AsRaster() alignment

[pracine]

I have a test raster that I display in Openjump:

CREATE OR REPLACE FUNCTION ST_TestRaster(h integer, w integer, val float8) 
    RETURNS raster AS 
    $$
    DECLARE
    BEGIN
        RETURN ST_AddBand(ST_MakeEmptyRaster(h, w, 0, 0, 1, 1, 0, 0, 0), '32BF', val, -1);
    END;
    $$
    LANGUAGE 'plpgsql';

SELECT ST_AsBinary((ST_PixelAsPolygons(ST_TestRaster(10, 10, 2))).geom);

Then I have a simple overlapping geometry:

SELECT ST_AsBinary(ST_Buffer(ST_MakePoint(8, 5), 4));

If I convert it to a raster aligned on the first raster I get:

SELECT ST_AsBinary((gv).geom), (gv).val 
FROM ST_PixelAsPolygons(ST_AsRaster(ST_Buffer(ST_MakePoint(8, 5), 4), 
                                    ST_TestRaster(10, 10, 2), 
                                    ST_BandPixelType(ST_TestRaster(10, 10, 2), 1), 
                                    1, 
                                    -10
                                   )
                        ) gv

The two rasters are well aligned. However if I try:

SELECT ST_MapAlgebraExpr(ST_TestRaster(10, 10, 2), 1, 
                         ST_AsRaster(ST_Buffer(ST_MakePoint(8, 5), 4), 
                                     ST_TestRaster(10, 10, 2), 
                                     ST_BandPixelType(ST_TestRaster(10, 10, 2), 1), 
                                     1, 
                                     -10
                                    ), 1, 
                         'rast1')

I get "NOTICE: The two rasters provided do not have the same alignment. Returning NULL"…

If I just replace the rasterized geometry with the original raster like this:

SELECT ST_MapAlgebraExpr(ST_TestRaster(10, 10, 2), 1, 
                         ST_TestRaster(10, 10, 2), 1, 
                         'rast1')

there is no problem…

I am trying to implement ST_Clip()…

[pracine]

Actually just this simple test fails:

SELECT ST_SameAlignment(ST_TestRaster(10, 10, 2), 
                        ST_AsRaster(ST_Buffer(ST_MakePoint(8, 5), 4), ST_TestRaster(10, 10, 2), '8BUI', 1, -10))

[Bborie Park]

If you check the individual metadata of each raster, you'll see that the Y-scale values differ. 1 vs -1. I'll add a check to ST_SameAlignment to compare the absolute values for scale as there is no difference between 1 and -1 in terms of scale, just direction.

[pracine]

Isn't the problem more in ST_AsRaster? It should just copy the exact same parameters from the raster.

[Bborie Park]

Yes and no. We make use of the GDAL's suggested extent, which likes a positive x-scale and and negative y-scale. I'll see if I can play with it.

[bnordgren]

Checking consistency of alignment may well be something that the "physically relevant geotransform" code/algorithm could handle #1331 . (pixel sizes are always positive). You just need to make a determination as to whether "same alignment" includes the case where the j basis vectors are opposed by 180 degrees. For instance, is it more important to you that the transform be very simple like a shift and add, or is it just necessary to have the grid cells be co-located?

In addition, it also becomes easier to circumscribe (not inscribe) gdal's suggested extent with a raster definition if the physically significant parameters are available. Circumscribing an extent (GBOX) with the grid definition of a source raster is exactly what motivated #1331. I have an algorithm to circumscribe a gbox with a raster worked out on paper.

This is actually the last known loose end for raster iterator.

[Bborie Park]

Can you try r8263, Pierre? I did some tweaking to the gdal rasterize code.

[pracine]

Fixed!