1

I want to calculate the proportion of each pixel value in a raster after clipping a much larger raster to my area of interest. However, after clipping and saving the raster from PostgreSQL database/PostGIS, the table created does not calculate the pixel number or proportion as expected.

Step 1: clip larger raster (Minnesota gridded soils) based on extent of "plotbounds" (polygon geometry table):

CREATE TABLE out.eonr_gaylord_2012_soils AS
SELECT ST_Clip(rast, p.geom)
FROM gssurgo_mn.mapunitraster_mn_10m, temp.eonr_gaylord_2012_plotbounds_data AS p 
WHERE ST_Intersects(p.geom, rast);

Step 2: Calculate the proportion of each pixel

SELECT DISTINCT pvc.VALUE, SUM(pvc.COUNT) AS tot_pix, 100 * SUM(pvc.COUNT) / SUM(SUM(pvc.COUNT)) OVER () AS prop
FROM out.eonr_gaylord_2012_soils, ST_ValueCount(rast, 1) AS pvc
GROUP BY pvc.VALUE, pvc.COUNT
ORDER BY pvc.VALUE

Output: Proportion of pixels incorrect

Why are there multiple rows (6) for each value?

I can load this raster into QGIS and it displays fine. I tried exporting to my local hard drive as a geotiff (from QGIS), then reloaded into PostgreSQL using raster2pgsql in the command prompt, and ran the same query to calculate pixel number and proportion (with the imported geotiff raster), and it calculates as expected:

Proportion of pixels correct

When looking at the Info from QGIS DB Manager, the original raster (the one that doesn't behave as expected) has a few differences from the imported raster (the one that works):

  1. Pages: 3 vs 1
  2. Rows: 128 vs 1
  3. Spatial ref: [none] vs "WGS 84 / UTM zone 15N (32615)"
  4. Constraints: [none] vs many (e.g., enforce_max_extent_rast, enforce_out_db_rast, etc.)

Version info:

SELECT version()

PostgreSQL 9.5.10, compiled by Visual C++ build 1800, 64-bit

SELECT PostGIS_full_version()

POSTGIS="2.3.5 r16110" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.2, released 2017/09/15" LIBXML="2.7.8" LIBJSON="0.12" RASTER

  • Is your raster tiled? In that case, you'd get a value per tile unless you used a union statement on the tiles, e.g. ST_ValueCount(ST_Union(rast), 1) – John Apr 12 at 5:11
  • @John Yes, it is. The union statement ends up working for me, although I do get an error saying "aggregate functions are not allowed in functions in FROM" if I use ST_ValueCount(ST_Union(rast), 1) in the FROM statement. This hint helped me find the solution to my problem though - thanks! – Tyler N Apr 15 at 20:00
0

Because my raster was tiled, it resulted in a calculation for each individual tile (thanks for the hint @John). Therefore, I can modify the query in Step 2 to achieve the desired result:

SELECT (pvc.valcount).VALUE, (pvc.valcount).COUNT AS tot_pix, (100.0 * (pvc.valcount).COUNT / SUM((pvc.valcount).COUNT) OVER ()) AS prop
FROM(
    SELECT (ST_ValueCount(rast_u.rast, 1)) AS valcount
    FROM (
        SELECT ST_Union(rast) AS rast
        FROM out.eonr_gaylord_2012_soils) AS rast_u) pvc

postgis union query

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.