[postgis-users] Masking raster values

[Pierre Racine]

Guido,

"no masks (you can create a mask as a band)" means that a mask is nothing different than a normal band. It's just a band that is mostly used to intersect (mask) another band having normal values. You can do that with... ... ST_Intersection(raster, raster)

So first you want to ST_MapAlgebra() your BQA to create a raster with 0 when those bits are set and 1 when they are not. You could use a SQL expression like "CASE [rast] & 61440 THEN 0 ELSE 1 END" in your ST_MapAlgebra() call.

CREATE TABLE mask AS
SELECT ST_MapAlgebraExpr(rast, bandnumber, "1BUI", "CASE [rast] & 61440 THEN 0 ELSE 1 END")
FROM yourLandsatTable

You can make this query faster by implementing a small callback PL/PGSQL function for the ST_MapAlgebra() variant taking a callback (http://postgis.net/docs/manual-2.1/RT_ST_MapAlgebra.html).

Second you want to intersect this "mask" with you some of your data band. Look at ST_Intersection(raster, raster).

CREATE TABLE newLandsat AS
SELECT ST_Intersection(a.rast, BAND_NUM, b.rast, 1, 'BAND1') rast
FROM yourLandsatTable a, mask b

You can do the tow operations in a single SQL statement but it is generally wiser and faster to make it in two.

As your question about creating raster operators I always found this idea very nice in theory but in practice you generally want to/have to be able to set a variety of options when doing operations involving many rasters: What should be the pixel type of the result? What to do with nodata values? Do you want the result to meet the extent of the first or the second raster, the union of them or the intersection of them? And so on... For sure you can establish default behavior for all of these in order to be able to construct nice and clean raster expressions, but then you become very restricted. PostGIS took the side of providing all the flexibility possible at the cost of sometimes hard to write function calls.

Hope this answer your questions.

Pierre

> -----Original Message-----
> From: postgis-users-bounces at lists.osgeo.org [mailto:postgis-users-
> bounces at lists.osgeo.org] On Behalf Of Guido LEMOINE
> Sent: Wednesday, February 12, 2014 12:19 PM
> To: postgis-users at lists.osgeo.org
> Subject: [postgis-users] Masking raster values
> 
> Dear List,
> 
> 
> 
> I am dabbling with some Landsat-8 imagery in PostGIS, using raster
> functionality. I am primarily interested in deriving polygon-delineated
> extracts from the various spectral bands (e.g. for a forest patch, an
> agricultural field, etc).
> 
> 
> 
> Landsat-8 provides the so-called BQA band (see
> http://landsat.usgs.gov/L8QualityAssessmentBand.php) which is a bit-
> coded set of quality parameters, for which the cloud and haze indicators
> (bits 12-15 (right to left)) are the most important (for me). I would want to
> mask out all pixels, in the spectral band(s), for which these bits are set in
> the BQA band.
> 
> 
> 
> Did anyone come across a solution that would address this kind of masking
> operation? The only reference I seem to be able to find is the statement
> that "no masks (you can create a mask as a band)" in the WKTRaster wiki
> page on osgeo, which is somewhat cryptic. I would know how to do this
> outside the database (using JAI), but it seems that masking is some core
> operation one would expect in raster functionality.
> 
> 
> 
> As a side issue, would it be possible (is it foreseen?) to create raster
> operators that work in an equivalent arithmetic way as on single variable
> (i.e. 2*rast, rast ~ 128, rast < 3) but produce rast as output (maybe through
> a mapping to ST_MapAlgebra functions)?
> 
> 
> 
> Any pointers welcome and excuses if I have overlooked discussions on this
> topic.
> 
> 
> 
> Guido Lemoine
> 
> 
> 
> Scientific Officer
> 
> European Commission, Joint Research Centre (JRC)
> 
> Institute for the Protection and Security of the Citizen (IPSC)
> 
> Global Security and Crisis Management Unit
> 
> 
> 
> Via E.Fermi 2749, I-21027 Ispra (VA) Italy, TP 268
> 
> Tel. +39 0332 786239 (direct line) Fax +39 0332 785154
> 
> e-mail: guido.lemoine at jrc.ec.europa.eu
> <mailto:guido.lemoine at jrc.ec.europa.eu>
> 
> web: http://globesec.jrc.ec.europa.eu <http://globesec.jrc.ec.europa.eu/>
> 
> 
> 
>