[raster] Real extent feature lost after metadata as views

[mloskot]

It seems, we've overlooked one feature while migrating from table to view for the raster_columns (#1216, #1319). It is quite an important feature.

With the raster_columns as table, it was possible to set real spatial extent of raster data, regardless of physical size of raster blob. It means, it was possible to 'cut off' any possible padding framing the real raster data.

With the raster_columns as view, spatial extent is calculated based on 'worldfile' metadata and actual width and height of raster matrix.

This is an undesired difference, generally, but seems to mostly affect uses of regular blocking.

Previously, it was possible to specify raster2pgsql with arbitrary block size. If the block size was not divinding the input raster evenly, the loader would generate padding along the right and bottom edges of right-most and bottom-most tiles, but the spatial extent of the whole tiles coverage written into raster_columns did not include that padding. It was extent of real data, ragardless physical pixel dimensions were increased by the padding.

Currently, with the raster_columns being a view, padding affects spatial extent.

For example:

• Input raster: ​utm.tif from GDAL test data:

f:\data\Test\pgraster\blocking>gdalinfo utm.tif
Driver: GTiff/GeoTIFF
Files: utm.tif
Size is 512, 512
Coordinate System is:
PROJCS["NAD27 / UTM zone 11N",
    GEOGCS["NAD27",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke 1866",6378206.4,294.9786982139006,
                AUTHORITY["EPSG","7008"]],
            AUTHORITY["EPSG","6267"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4267"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26711"]]
Origin = (440720.000000000000000,3751320.000000000000000)
Pixel Size = (60.000000000000000,-60.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  440720.000, 3751320.000) (117d38'28.21"W, 33d54' 8.47"N)
Lower Left  (  440720.000, 3720600.000) (117d38'20.79"W, 33d37'31.04"N)
Upper Right (  471440.000, 3751320.000) (117d18'32.07"W, 33d54'13.08"N)
Lower Right (  471440.000, 3720600.000) (117d18'28.50"W, 33d37'35.61"N)
Center      (  456080.000, 3735960.000) (117d28'27.39"W, 33d45'52.46"N)
Band 1 Block=512x16 Type=Byte, ColorInterp=Gray

• Block 256x256

raster2pgsql -C -r -t 256x256 public.utm256 utm.tif

The loader cut the raster evenly, so no padding added, spatial extent the same as original:

test=# SELECT ST_Extent(rast::geometry) FROM utm256;
             st_extent
------------------------------------
 BOX(440720 3720600,471440 3751320)
(1 row)

• Block 300x300

raster2pgsql -C -r -t 300x300 public.utm300 utm.tif

The loader added padding, so spatial extent increased rightwards and downwards:

test=# SELECT ST_Extent(rast::geometry) FROM utm300x300;
             st_extent
------------------------------------
 BOX(440720 3715320,476720 3751320)
(1 row)

This is a significant change in behaviour, a regression which I'm reporting as a bug.

Apparently, PostGIS Raster driver in GDAL 1.9.2 is affected. I tried to export the coverages from the two tables loaded above, on Windows:

set CONN=PG:dbname=test user=mloskot mode=2
FOR %%t IN (utm256 utm300) DO (
    gdal_translate -of GTIff "%CONN% table=%%t" %%t.tif
	gdalinfo %%t.tif
)

First, similarly to the SQL results above, gdalinfo reports changed spatial extent for the utm300.tif file. Second, the padding strips are exported, whereas in normal situation they should be discarded. I'm going to attach both files, utm256.tif and utm300.tif:

---

I'd like to bring the previous behaviour back.

I'd like to propose change in the header of serialised format, and add 4 slots of 16 bytes:

struct rt_raster_serialized_t {
    /*---[ 8 byte boundary ]---{ */
    uint32_t size; /* required by postgresql: 4 bytes */
    uint16_t version; /* format version (this is version 0): 2 bytes */
    uint16_t numBands; /* Number of bands: 2 bytes */

    /* }---[ 8 byte boundary ]---{ */
    double scaleX; /* pixel width: 8 bytes */

    /* }---[ 8 byte boundary ]---{ */
    double scaleY; /* pixel height: 8 bytes */

    /* }---[ 8 byte boundary ]---{ */
    double ipX; /* insertion point X: 8 bytes */

    /* }---[ 8 byte boundary ]---{ */
    double ipY; /* insertion point Y: 8 bytes */

    /* }---[ 8 byte boundary ]---{ */
    double skewX; /* skew about the X axis: 8 bytes */

    /* }---[ 8 byte boundary ]---{ */
    double skewY; /* skew about the Y axis: 8 bytes */

    /* }---[ 8 byte boundary ]--- */
    int32_t srid; /* Spatial reference id: 4 bytes */
    uint16_t width; /* pixel columns: 2 bytes */
    uint16_t height; /* pixel rows: 2 bytes */
	
	/* }---[ 8 byte boundary ]--- */
    uint16_t upperLeftX; /* first column of raster data, Zero if no padding leftwards */
    uint16_t upperLeftY; /* first row of raster data, Zero if no padding upwards */
    uint16_t lowerRightX; /* last column of raster data, width if no padding rightwards */
    uint16_t lowerRightY; /* last row of raster data, height if no padding downwards */
};

I'm going to post this issue to postgis-devel too and request for comments.

[mloskot]

Loaded using raster2pgsql -C -r -t 256x256 public.utm256 utm.tif, then exported with gdal_translate to single file

[mloskot]

Loaded using raster2pgsql -C -r -t 300x300 public.utm300 utm.tif, then exported with gdal_translate to single file

[mloskot]

The ​Real extent feature lost after metadata as views thread consists of detailed discussion of the problem.

[mloskot]

Here is the final proposal for the updated regular blocking specification proposed by Pierre after the course of lengthy discussion: ​Regular raster blocking

Here is ​confirmation of agreement with minor editorial suggestions. Namely, the three rules are named as follows:

1. Tile raster metadata
2. Tile raster location
3. Tile raster dimensions

So, we have the final spec of the regular blocking. Further editing and polishing should happen in the SVN.

Pierre, thanks for your support!

[mloskot]

An intermediate correction to blocksize constraints that replaces '=' with '<=' operator for tile dimension and blocksize comparison. NOTE: this is not complete implementation of the rule 3) from Pierr's proposal because it does not restrict '<=' comaprisons to the right-most and bottom tiles only. It's an intermediate enabler to allow the use of PostGIS Raster sooner.

[Bborie Park]

I'm fine with the proposal "on paper". Implementing is a whole different nightmare as parts of it may not be able to be implemented, specifically testing that any two tiles can only touch.

[pracine]

Replying to dustymugs:

I'm fine with the proposal "on paper". Implementing is a whole different nightmare as parts of it may not be able to be implemented, specifically testing that any two tiles can only touch.

We removed the need to test for gaps and overlaps. Everything is now just based on metadatas, upper left corners and tile sizes. I don't think this involves relationships like touch or intersects anymore.

I think the test of alignment of upper left corners resume to something very close to ST_SameAlignment(). I am more worry about how to determine what tiles are the right-most and bottom ones.

[Bborie Park]

Closing ticket as the mess here is broken down into and resolved through individual tickets...

#2148 #2149 #2150 #2069 #826 #2143