[postgis-users] Masking raster values

Fri Feb 14 00:59:13 PST 2014

Hey,
short suggestion,
you may want to try to not use a "case" for perf reasons, in perticular if
it is callen often
You could do something equivalent like this :

floor(rast/61440).
If rast is an int you may even skip the floor (don't have postgres here to
test it).

Cheers,

Rémi-C

2014-02-13 21:42 GMT+01:00 guido lemoine <guido.lemoine at jrc.ec.europa.eu>:

> Pierre,
>
> Thanks very much for the answer. Your suggestion works with some
> modification.
> I had to modify to
>
> CASE WHEN [rast] > 61440 THEN 0 ELSE 1 END
>
> for the Expr in MapAlgebraExpr.
>
> (the WHEN was missing, and the logical and (&) does not work).
>
> I use tiled raster storage with 64 by 64 tile size, and the single query
> version combining with ST_Intersection works fast enough.
> The bands are from the same image, so all nicely co-located (I am storing
> bands separately, but may reconsider that).
>
> For the operator issue, I will dig a little deeper.
>
> Guido
>
>
>
> On 02/13/14, postgis-users-request at lists.osgeo.org wrote:
>
> Send postgis-users mailing list submissions to
> postgis-users at lists.osgeo.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
> or, via email, send a message with subject or body 'help' to
> postgis-users-request at lists.osgeo.org
>
> You can reach the person managing the list at
> postgis-users-owner at lists.osgeo.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of postgis-users digest..."
>
> ------------------------------
>
> Today's Topics:
>
>    1. Re: 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-<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 <guido.lemoine at jrc.ec.europa.eu>>
>
> >
> > web: http://globesec.jrc.ec.europa.eu <http://globesec.jrc.ec.europa.eu/
> >
> >
> >
> >
> >
>
>
> ------------------------------
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20140214/1d3ee4ea/attachment.html>