1

Does anyone know a nice and easy way to recenter a tiled raster table in PostGIS?

By recentering I mean that I have a raster that is centered on the Pacific, but I need to recenter it on Greenwhich/Europe. In a sense the effect would be to offset pixels while "rolling" pixels that fall off the right edge back around to the left edge. It would be easier if the entire raster data was just a single big raster, but due to large memory size I have chopped up the raster into several subrasters/tiles with the raster2pgsql -t option (which is the PostGIS approach to raster data), which means the operation would have to be done on a tile by tile basis making the solution slightly more challenging.

1 Answer 1

1

I eventually developed a solution and put it into a reusable function. It is pretty flexible: regardless of the coordinate bounds or system of the raster, you simply offset it rightwards by how many coordinates you wish, and all whole and partial tiles that fall off the right side will be clipped and wrapped around to the left side.

-- Function for horizontally rolling a tiled raster table by x coordinates rightwards and wrapping the parts that fall off back around to the left side
-- MAYBE TODO: Make faster by instead of doing clip on all tiles, do a "UNION ALL" type combine on:
-- 1) those tiles not intersecting the new xcenter vertical line, 
-- 2) those tiles intersecting clipped on the left side, and 
-- 3) those tiles intersecting clipped on the right side.

CREATE FUNCTION Roll_Raster_Table(text, text, numeric) 
RETURNS TABLE(rast raster) AS  
$$
DECLARE
    _rastertablename ALIAS FOR $1;
    _rastercolname ALIAS FOR $2;
    _xoffset ALIAS FOR $3;
BEGIN
    RETURN QUERY EXECUTE '
    WITH
    _wholerastmeta AS
        -- whole raster table meta (assumes right and up increasing orientation of coords)
        (
        WITH 
            leftedges AS (SELECT ST_UpperLeftX('||_rastercolname||') AS val FROM '||_rastertablename||'),
            rightedges AS (SELECT ST_UpperLeftX('||_rastercolname||')+ST_Width('||_rastercolname||')*ST_ScaleX('||_rastercolname||') AS val FROM '||_rastertablename||'),
            topedges AS (SELECT ST_UpperLeftY('||_rastercolname||') AS val FROM '||_rastertablename||'),
            bottomedges AS (SELECT ST_UpperLeftY('||_rastercolname||')+ST_Height('||_rastercolname||')*ST_ScaleY('||_rastercolname||') AS val FROM '||_rastertablename||')
        SELECT MIN(leftedges.val) AS leftedge, MAX(rightedges.val) AS rightedge, MIN(bottomedges.val) AS bottomedge, MAX(topedges.val) AS topedge
        FROM leftedges, rightedges, bottomedges, topedges
        )
        ,
    _clipped AS
        -- clip all tiles to the left and right of new rightedge
        (
        SELECT ST_Clip('||_rastercolname||', ST_SetSRID(ST_MakeEnvelope(leftedge, bottomedge, rightedge-'||_xoffset||', topedge), 4326)) AS rast
        FROM '||_rastertablename||', _wholerastmeta
        WHERE '||_rastercolname||' && ST_MakeEnvelope(leftedge, bottomedge, rightedge-'||_xoffset||', topedge)

        UNION ALL

        SELECT ST_Clip('||_rastercolname||', ST_SetSRID(ST_MakeEnvelope(rightedge-'||_xoffset||', bottomedge, rightedge, topedge), 4326)) AS rast
        FROM '||_rastertablename||', _wholerastmeta
        WHERE '||_rastercolname||' && ST_MakeEnvelope(rightedge-'||_xoffset||', bottomedge, rightedge, topedge)
        )
    -- left offset any tiles that are beyond the new rightedge, and rightoffset the ones on the left
    SELECT (CASE WHEN ST_UpperLeftX(rast) >= rightedge-'||_xoffset||'
        THEN ST_SetUpperLeft(rast, ST_UpperLeftX(rast)-(rightedge-leftedge)+'||_xoffset||', ST_UpperLeftY(rast))
        ELSE ST_SetUpperLeft(rast, ST_UpperLeftX(rast)+'||_xoffset||', ST_UpperLeftY(rast))
        END) AS rast
    FROM _clipped, _wholerastmeta
    ';
END
$$ language 'plpgsql';  

Note however, for some reason it seems to run slightly slow in some cases, so if you need something faster and all you need is to recenter a global pacific-centered raster to Europe/Africa, then here is a hardcoded query you can use:

WITH 
_input AS
    (SELECT rast FROM rastertable)
    ,
_clipped AS
-- clip all tiles to the left and right of new rightedge
    (       
    SELECT ST_Clip(rast, ST_SetSRID(ST_MakeEnvelope(0, -90, 180, 90), 4326)) AS rast
    FROM _input
    WHERE rast && ST_MakeEnvelope(0, -90, 180, 90) --intersects left side

    UNION ALL

    SELECT ST_Clip(rast, ST_SetSRID(ST_MakeEnvelope(180, -90, 360, 90), 4326)) AS rast
    FROM _input
    WHERE rast && ST_MakeEnvelope(180.001, -90, 360, 90) --intersects right side (180.000001 to avoid weird edge intersection empty band clip error)
    )
-- switch sides the tiles around the middle (left offset any tiles that are beyond the middle, and right offset all the other ones)
SELECT (CASE WHEN ST_UpperLeftX(rast) >= 180
    THEN ST_SetUpperLeft(rast, ST_UpperLeftX(rast)-180, ST_UpperLeftY(rast))
    ELSE ST_SetUpperLeft(rast, ST_UpperLeftX(rast)+180, ST_UpperLeftY(rast))
    END) AS rast
FROM _clipped

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.