<div dir="ltr"><div><div><div>Hey,<br>short suggestion, <br>you may want to try to not use a "case" for perf reasons, in perticular if it is callen often<br></div>You could do something equivalent like this :<br>
<br></div>floor(rast/61440).<br></div>If rast is an int you may even skip the floor (don't have postgres here to test it).<br><br>Cheers,<br><br>Rémi-C<br></div><div class="gmail_extra"><br><br><div class="gmail_quote">
2014-02-13 21:42 GMT+01:00 guido lemoine <span dir="ltr"><<a href="mailto:guido.lemoine@jrc.ec.europa.eu" target="_blank">guido.lemoine@jrc.ec.europa.eu</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div style="font-family:'Times New Roman';font-size:16px"><div>Pierre,</div><div><br></div><div>Thanks very much for the answer. Your suggestion works with some modification. </div><div>I had to modify to </div><div>
<br></div><div>CASE WHEN [rast] > 61440 THEN 0 ELSE 1 END</div><div><br></div><div>for the Expr in MapAlgebraExpr.</div><div><br></div><div>(the WHEN was missing, and the logical and (&) does not work).</div><div><br>
</div><div>I use tiled raster storage with 64 by 64 tile size, and the single query version combining with ST_Intersection works fast enough.</div><div>The bands are from the same image, so all nicely co-located (I am storing bands separately, but may reconsider that).</div>
<div><br></div><div>For the operator issue, I will dig a little deeper.</div><div><br></div><div>Guido</div><div><br></div><div><br></div><div><br></div><div>On 02/13/14, <a href="mailto:postgis-users-request@lists.osgeo.org" target="_blank">postgis-users-request@lists.osgeo.org</a> wrote:</div>
<blockquote style="border-left:1px solid #00f;padding-left:13px;margin-left:0" type="cite"><div>Send postgis-users mailing list submissions to<br>        <a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><br>
<br>To subscribe or unsubscribe via the World Wide Web, visit<br>   <a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br>
or, via email, send a message with subject or body 'help' to<br>  <a href="mailto:postgis-users-request@lists.osgeo.org" target="_blank">postgis-users-request@lists.osgeo.org</a><br><br>You can reach the person managing the list at<br>
        <a href="mailto:postgis-users-owner@lists.osgeo.org" target="_blank">postgis-users-owner@lists.osgeo.org</a><br><br>When replying, please edit your Subject line so it is more specific<br>than "Re: Contents of postgis-users digest..."<br>
<p></p><hr size="2"><p>Today's Topics:<br><br>   1. Re: Masking raster values (Pierre Racine)<br></p><hr size="2"><p></p><div class="">Guido,<br><br>"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)<br>
<br>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.<br>
<br>CREATE TABLE mask AS<br>SELECT ST_MapAlgebraExpr(rast, bandnumber, "1BUI", "CASE [rast] & 61440 THEN 0 ELSE 1 END")<br>FROM yourLandsatTable<br><br>You can make this query faster by implementing a small callback PL/PGSQL function for the ST_MapAlgebra() variant taking a callback (<a href="http://postgis.net/docs/manual-2.1/RT_ST_MapAlgebra.html)." target="_blank">http://postgis.net/docs/manual-2.1/RT_ST_MapAlgebra.html).</a><br>
<br>Second you want to intersect this "mask" with you some of your data band. Look at ST_Intersection(raster, raster).<br><br>CREATE TABLE newLandsat AS<br>SELECT ST_Intersection(a.rast, BAND_NUM, b.rast, 1, 'BAND1') rast<br>
FROM yourLandsatTable a, mask b<br><br>You can do the tow operations in a single SQL statement but it is generally wiser and faster to make it in two.<br><br>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.<br>
<br>Hope this answer your questions.<br><br>Pierre<br><br>> -----Original Message-----<br></div><div><div class="h5">> From: <a href="mailto:postgis-users-bounces@lists.osgeo.org" target="_blank">postgis-users-bounces@lists.osgeo.org</a> [mailto:<a href="mailto:postgis-users-" target="_blank">postgis-users-</a> <postgis-users-><br>
> <a href="mailto:bounces@lists.osgeo.org" target="_blank">bounces@lists.osgeo.org</a>] On Behalf Of Guido LEMOINE<br>> Sent: Wednesday, February 12, 2014 12:19 PM<br>> To: <a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><br>
> Subject: [postgis-users] Masking raster values<br>> <br>> Dear List,<br>> <br>> <br>> <br>> I am dabbling with some Landsat-8 imagery in PostGIS, using raster<br>> functionality. I am primarily interested in deriving polygon-delineated<br>
> extracts from the various spectral bands (e.g. for a forest patch, an<br>> agricultural field, etc).<br>> <br>> <br>> <br>> Landsat-8 provides the so-called BQA band (see<br>> <a href="http://landsat.usgs.gov/L8QualityAssessmentBand.php)" target="_blank">http://landsat.usgs.gov/L8QualityAssessmentBand.php)</a> which is a bit-<br>
> coded set of quality parameters, for which the cloud and haze indicators<br>> (bits 12-15 (right to left)) are the most important (for me). I would want to<br>> mask out all pixels, in the spectral band(s), for which these bits are set in<br>
> the BQA band.<br>> <br>> <br>> <br>> Did anyone come across a solution that would address this kind of masking<br>> operation? The only reference I seem to be able to find is the statement<br>> that "no masks (you can create a mask as a band)" in the WKTRaster wiki<br>
> page on osgeo, which is somewhat cryptic. I would know how to do this<br>> outside the database (using JAI), but it seems that masking is some core<br>> operation one would expect in raster functionality.<br>> <br>
> <br>> <br>> As a side issue, would it be possible (is it foreseen?) to create raster<br>> operators that work in an equivalent arithmetic way as on single variable<br>> (i.e. 2*rast, rast ~ 128, rast < 3) but produce rast as output (maybe through<br>
> a mapping to ST_MapAlgebra functions)?<br>> <br>> <br>> <br>> Any pointers welcome and excuses if I have overlooked discussions on this<br>> topic.<br>> <br>> <br>> <br>> Guido Lemoine<br>> <br>
> <br>> <br>> Scientific Officer<br>> <br>> European Commission, Joint Research Centre (JRC)<br>> <br>> Institute for the Protection and Security of the Citizen (IPSC)<br>> <br>> Global Security and Crisis Management Unit<br>
> <br>> <br>> <br>> Via E.Fermi 2749, I-21027 Ispra (VA) Italy, TP 268<br>> <br>> Tel. <a href="tel:%2B39%200332%20786239" value="+390332786239" target="_blank">+39 0332 786239</a> (direct line) Fax <a href="tel:%2B39%200332%20785154" value="+390332785154" target="_blank">+39 0332 785154</a><br>
> <br>> e-mail: <a href="mailto:guido.lemoine@jrc.ec.europa.eu" target="_blank">guido.lemoine@jrc.ec.europa.eu</a><br></div></div>> <mailto:<a href="mailto:guido.lemoine@jrc.ec.europa.eu" target="_blank">guido.lemoine@jrc.ec.europa.eu</a> <<a href="mailto:guido.lemoine@jrc.ec.europa.eu" target="_blank">guido.lemoine@jrc.ec.europa.eu</a>>><div class="">
<br>> <br>> web: <a href="http://globesec.jrc.ec.europa.eu" target="_blank">http://globesec.jrc.ec.europa.eu</a> <<a href="http://globesec.jrc.ec.europa.eu/" target="_blank">http://globesec.jrc.ec.europa.eu/</a>><br>
> <br>> <br>> <br>> <br><br><br></div><p></p><div class=""><hr size="2"><p>_______________________________________________<br>postgis-users mailing list<br><a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a></p></div></div></blockquote></div>
<br>_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br></blockquote></div><br></div>