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, [email protected] wrote: > > Send postgis-users mailing list submissions to > [email protected] > > 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 > [email protected] > > You can reach the person managing the list at > [email protected] > > 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: [email protected] [mailto:postgis-users- > > <postgis-users-> > > [email protected]] On Behalf Of Guido LEMOINE > > Sent: Wednesday, February 12, 2014 12:19 PM > > To: [email protected] > > 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: [email protected] > > <mailto:[email protected] <[email protected]>> > > > > web: http://globesec.jrc.ec.europa.eu <http://globesec.jrc.ec.europa.eu/> > > > > > > > > > > > > > _______________________________________________ > postgis-users mailing list > [email protected] > http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users >
_______________________________________________ postgis-users mailing list [email protected] http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
