ST_SummaryStats gives unknown values (infinity and different values) when clipping raster with buffer

Question

I am loading a raster to postgis (see postgis version and command below)

raster2pgsql /path_to_raster/temp-221025_12.tif -d icon_test -F -n name | PGPASSWORD=<PASS> psql -d <myDB> -h <host> -U <user> -p <port>

POSTGIS="3.1.2 cbe925d" [EXTENSION] PGSQL="120" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.1.2, released 2016/10/24" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.2.1" WAGYU="0.5.0 (Internal)" RASTER

I am trying to get stats of pixels near a point with this query:

SELECT name,ST_SummaryStats(
ST_Clip(rast, ST_buffer(ST_GeomFromText('POINT(-90.37083 15.47083)',4326), 3000*0.00001)) 
) from icon_test

But it gives me BAD stats according to my knowledge (check columns in the image)

I have loaded the tif to QGIS and overlaying the point buffer and it overlays to 4 pixels as the stats says but I was expecting mean as ((18.797+18.391+17.976+17.544)/4) = 18.177 and max=18.797 and min=17.544 (based on the image)

Why ST_SummaryStats is not giving the expected stats ?

Link to raster if you want to download

UPDATE 1 using envelope works for this coord:

SELECT name,
(ST_SummaryStats(
        ST_Clip(rast, st_envelope(ST_buffer(ST_GeomFromText('POINT(-90.37083 15.47083)',4326), 3000*0.00001))
        ) 
    )
).* from icon_test

Result:

name              |count|sum              |mean             |stddev            |min               |max               |
------------------+-----+-----------------+-----------------+------------------+------------------+------------------+
temp-221025_12.tif|    4|72.70815658569336|18.17703914642334|0.4666503646805367|17.544031143188477|18.797155380249023|

Result for coord. POINT(-90.6 14.52):

name              |count|sum      |mean     |stddev  |min                                                                                                                                                                                                                                                            |max               |
------------------+-----+---------+---------+--------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+------------------+
temp-221025_12.tif|    4|-Infinity|-Infinity|Infinity|-17976931348623157000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000|13.930984497070312|

UPDATE 2 Raster STATS same with POSTGIS and cli-GDAL:

count |sum             |mean              |stddev            |min               |max               |
------+----------------+------------------+------------------+------------------+------------------+
990525|23331638.4433757|23.554820366346835|5.6419397045642325|-4.649328708648682|31.385828018188477|

UPDATE 3 pixel values of clip:

select 
 ST_PixelAsPoints(
    ST_CLIP(rast, ST_buffer(ST_GeomFromText('POINT (-90.37083 15.47083)',4326), 3000*0.00001))
    ) geom
from icon_test

geom                                                                             |
---------------------------------------------------------------------------------+
(0101000020E610000033333333339B56C0CCCCCCCCCC0C2F40,-1.7976931348623157e+308,1,1)|
(0101000020E610000000000000009856C0CCCCCCCCCC0C2F40,-1.7976931348623157e+308,2,1)|
(0101000020E610000033333333339B56C03233333333F32E40,-1.7976931348623157e+308,1,2)|
(0101000020E610000000000000009856C03233333333F32E40,17.544031143188477,2,2)      |

UPDATE 4: I am unable to get this to work for all rasters, any other ideas?

Answer 1

I'm not 100% sure, but I think the way ST_Clip works is that it keep the data if and only if the center of the pixel is in your geom.

So seeing your image only the bottom right pixel will be kept. But the resulting raster would still include 4 pixels because the geom intersects with 4 pixels. Others pixels will be with the value 'nodata', which by default will be ST_MinPossibleValue(ST_BandPixelType(rast, band)), in your case probably -infinity.

One way to handle that would probably be to add half of the pixel size to the buffer (if you have square pixels, otherwise it will be more complicated), but you need to handle the corners too. So I think you need to switch to square selection after having your circle:

SELECT 
  name,
  ST_SummaryStats(
    ST_Clip(
      rast, 
      ST_Enveloppe(
        ST_buffer(
          ST_GeomFromText('POINT(-90.37083 15.47083)', 4326),
          3000*0.00001 + ST_PixelWidth(rast)/2
        )
      )
    ) 
  ) 
from icon_test

It's just an idea, maybe there is specific cases that I missed.

And after that, you still need to handle missing values in your resulting raster. Maybe there is a problem with the nodata_value of the resulting raster, because it should be ignored by ST_SummaryStats but it doesn't seems like it is. Maybe try to set a nodataval in ST_Clip ?

EDIT:

Like I said, I think you second problem is the nodata value, you need to set it for it to work as intended. I think this should work:

SELECT 
  name,
  ST_SummaryStats(
    ST_SetBandNoDataValue(
      ST_Clip(
        rast, 
        ST_Enveloppe(
          ST_buffer(
            ST_GeomFromText('POINT(-90.37083 15.47083)', 4326),
            3000*0.00001 + ST_PixelWidth(rast)/2
          )
        ),
        FLOAT8 '-infinity'
      ),
      FLOAT8 '-infinity'
    ),
    TRUE
  ) 
from icon_test