1

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) enter image description here

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) enter image description here

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?

1
  • 1
    Like I said, the nodata value needs to be set, ohterwise you're doing stat on -infinity values. Can you try the edit I did on my answer ? Nov 2 at 10:21

1 Answer 1

0

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
5
  • tested with envelope and works for the testing point but not for other points, see UPDATE 1 of my question. Stats of raster as geotiff with GDAL and with Postgis shows same values. seems no missing or novalues points.
    – neavilag
    Oct 27 at 14:14
  • Instead of ST_SummaryStats, can you display the values of the raster returned by ST_Clip to see what's going on ? Oct 27 at 14:25
  • is what I attempted to do and got this... and seems that is the data take for summary but pixels has data if I enlarge buffer got more valid values but always the -1.79e+318 values... see UPDATE 3
    – neavilag
    Oct 27 at 14:57
  • You can surround ST_PixelAsPoints() with ().* for better display: (ST_PixelAsPoints()).* Oct 27 at 17:49
  • 2
    -1.7976931348623157e+308 is minus infinity, so the result of ST_SummaryStats is logical. Like I said, you get the pixel but filled with nodata, but because nodata is not set in your original raster you get -infinty. You need to make ST_SummaryStats understand that the -infinity is nodata. Try using ST_SetBandNoDataValue to set nodata to -infinity. You can also use your own nodata (like -10000 for example) and set it in both ST_SetBandNoDataValue and ST_Clip. Oct 28 at 7:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.