[raster] inconsistent nodata values returned from multiband raster

[gdrolet]

Inconsistent nodata values are returned when calling ST_Value on an out-db multiband raster:

SELECT  gid,
        (ST_WorldToRasterCoord(rast, geom)).columnx x,
	(ST_WorldToRasterCoord(rast, geom)).rowy y,
	ST_Value(rast, 1, geom, FALSE) sd1,
	ST_Value(rast, 2, geom, FALSE) sd2,
 	ST_Value(rast, 3, geom, FALSE) sd3,
	ST_Value(rast, 4, geom, FALSE) sd4,
	ST_Value(rast, 5, geom, FALSE) sd5,
	ST_Value(rast, 6, geom, FALSE) sd6	
FROM sg.raster_10x10_mult b
	INNER JOIN (
		SELECT  gid, geom
		FROM    sg.points
		limit	10
	) l
	ON ST_Intersects(rast, geom)
ORDER BY gid;

Pixels at locations with no data should all return -9999: bands 3 to 6 return 0.

SELECT  x band_number, 		
       (ST_BandMetadata(rast, x)).nodatavalue
FROM   sg.raster_10x10_mult,(
         SELECT generate_series(1, 6) x
       ) foo			
WHERE  rid = 1; 

 band_number | nodatavalue
-------------+-------------
           1 |       -9999
           2 |       -9999
           3 |       -9999
           4 |       -9999
           5 |       -9999
           6 |       -9999
(6 rows)

Raster sg.raster_10x10_mult was loaded using:

raster2pgsql -s 3175 -p -R -P -F -e -Y -t 10x10 BLD_all_M.tif sg.raster_10x10_mult | psql -U postgres -d testdb -h localhost

> gdalinfo BLD_all_M.tif 

Driver: GTiff/GeoTIFF
Files: BLD_all_M_1km.tif
Size is 2659, 3240
Coordinate System is:
PROJCS["NAD83 / Great Lakes and St Lawrence Albers",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",42.122774],
    PARAMETER["standard_parallel_2",49.01518],
    PARAMETER["latitude_of_center",45.568977],
    PARAMETER["longitude_of_center",-83.248627],
    PARAMETER["false_easting",1000000],
    PARAMETER["false_northing",1000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","3175"]]
Origin = (1188631.738954391100000,2982399.663283404900000)
Pixel Size = (801.504881290890810,-801.504881290890810)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 1188631.739, 2982399.663) ( 79d42'32.91"W, 63d42'27.63"N)
Lower Left  ( 1188631.739,  385523.848) ( 81d 2'40.85"W, 40d 0'56.80"N)
Upper Right ( 3319833.218, 2982399.663) ( 43d19'55.20"W, 58d11'52.36"N)
Lower Right ( 3319833.218,  385523.848) ( 57d 6' 1.42"W, 36d35'24.78"N)
Center      ( 2254232.479, 1683961.756) ( 65d28'29.73"W, 50d27'54.43"N)
Band 1 Block=2659x1 Type=Int16, ColorInterp=Gray
  NoData Value=-9999
Band 2 Block=2659x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-9999
Band 3 Block=2659x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-9999
Band 4 Block=2659x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-9999
Band 5 Block=2659x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-9999
Band 6 Block=2659x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-9999

PostGIS version:

POSTGIS="2.3.0dev r14604" GEOS="3.5.0-CAPI-1.9.0 r4090" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 2.0.1, released 2015/09/15" LIBXML="2.7.8" LIBJSON="0.12" RASTER

[gdrolet]

Here's the output from the point/raster overlay query (not included in my original post):

 gid | x | y  |  sd1  |  sd2  | sd3  | sd4  | sd5  | sd6
-----+---+----+-------+-------+------+------+------+------
   1 | 5 |  8 |  1316 |  1335 | 1366 | 1416 | 1477 | 1531
   2 | 9 | 10 |  1241 |  1259 | 1290 | 1341 | 1401 | 1455
   3 | 3 |  6 |  1230 |  1249 | 1280 | 1330 | 1391 | 1445
   4 | 1 |  8 | -9999 | -9999 |    0 |    0 |    0 |    0
   5 | 1 |  6 |  1218 |  1237 | 1267 | 1318 | 1379 | 1433
   7 | 1 |  3 |  1253 |  1271 | 1302 | 1353 | 1413 | 1467
   8 | 3 |  3 |  1351 |  1369 | 1400 | 1451 | 1511 | 1566
   9 | 8 |  2 | -9999 | -9999 |    0 |    0 |    0 |    0
  10 | 2 |  6 |  1188 |  1207 | 1237 | 1288 | 1349 | 1403
(9 rows)

[robe]

Do you have a sample raster file in question here. Given the results, there is not enough for me to tell if anything is wrong here? The results look perfectly fine to me.

In your first query they returned all -9999 as expected.

In your overlay, you are looking at a tile that intersects your geometry and are pulling out a value of the cell that intersects your point.

All that tells me is yes, at least one of the bands of each row should return a non-nodata value, which it does for bands 3-6 (value of 0) which is not nodata.

For now I'm going to close this out as an invalid, because I suspect nothing is wrong here.