[prev in list] [next in list] [prev in thread] [next in thread]
List: postgis-users
Subject: Re: [postgis-users] [Raster] Finding pixels intersecting an extent (polygon)
From: Jean Marchal <jean.d.marchal () gmail ! com>
Date: 2015-02-06 17:10:12
Message-ID: CACF6B3W2uXsOfKA88uG0kgvcPKdHhntOVbh1pMDaFfoQrRiLgQ () mail ! gmail ! com
[Download RAW message or body]
[Attachment #2 (multipart/alternative)]
For the posterity, I modified very slightly the query by changing the
returnband argument of ST_Intersection to return only one band.
WITH foo AS (
SELECT
ST_AsRaster(
ST_GeomFromText('POLYGON ((-52.54178994517749 46.99199259385565,
-52.54178994517749 46.996897959881, -52.53436080387823 46.996897959881,
-52.53436080387823 46.99199259385565, -52.54178994517749
46.99199259385565))', 4269),
rast,
'8BUI',
touched := True
) AS rast
FROM elev
LIMIT 1
)
SELECT
ST_Intersection(
elev.rast,
foo.rast,
'BAND1'
) AS rast
FROM elev
JOIN foo
ON ST_Intersects(elev.rast, foo.rast)
Thanks again,
Jean
2015-02-05 15:29 GMT-08:00 Jean Marchal <jean.d.marchal@gmail.com>:
> It worked! Thanks a lot!
>
> Jean
>
> 2015-02-05 15:18 GMT-08:00 Bborie Park <dustymugs@gmail.com>:
>
> Try something like:
>>
>> WITH foo AS (
>> SELECT
>> ST_AsRaster(
>> ST_GeomFromText('POLYGON ((-52.54178994517749 46.99199259385565,
>> -52.54178994517749 46.996897959881, -52.53436080387823 46.996897959881,
>> -52.53436080387823 46.99199259385565, -52.54178994517749
>> 46.99199259385565))', 4269),
>> rast,
>> '8BUI',
>> touched := True
>> ) AS rast
>> FROM elev
>> LIMIT 1
>> )
>> SELECT
>> ST_Intersection(
>> elev.rast,
>> foo.rast
>> ) AS rast
>> FROM elev
>> JOIN foo
>> ON ST_Intersects(elev.rast, foo.rast)
>>
>> The idea is to explicitly convert the geometry into a raster using one of
>> the elev rasters as a reference. The "touched := True" slightly changes the
>> behavior of the rasterization.
>>
>> From the docs:
>>
>> The optional touched parameter defaults to false and maps to the GDAL
>> ALL_TOUCHED rasterization option, which determines if pixels touched by
>> lines or polygons will be burned.
>>
>> On Thu, Feb 5, 2015 at 2:45 PM, Jean Marchal <jean.d.marchal@gmail.com>
>> wrote:
>>
>>> Bborie,
>>>
>>> I am running PostGIS 2.1.5. Here is the output of the
>>> postgis_full_version():
>>>
>>> "POSTGIS="2.1.5 r13152" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23
>>> September 2009" GDAL="GDAL 1.9.0, released 2011/12/29" LIBXML="2.8.0"
>>> LIBJSON="UNKNOWN" RASTER"
>>>
>>> My query looks like this:
>>>
>>> SELECT ST_Intersection(ST_GeomFromText('POLYGON ((-52.54178994517749
>>> 46.99199259385565, -52.54178994517749 46.996897959881, -52.53436080387823
>>> 46.996897959881, -52.53436080387823 46.99199259385565, -52.54178994517749
>>> 46.99199259385565))', 4269), rast) as rast
>>> FROM elev
>>>
>>> David,
>>>
>>> This query does not returns all the polygons that intersect my polygon.
>>> I think it returns only those where the centroid is inside the polygon or
>>> is it covered based? (like 50% of the pixel intersect with the polygon).
>>>
>>> Thanks for your rapid answers!
>>>
>>> Jean
>>>
>>> 2015-02-05 14:31 GMT-08:00 David Haynes <haynesd2@gmail.com>:
>>>
>>> select ST_Clip(r.rast,p.geom) as rast
>>>> from polygon p inner join raster r on ST_intersects(r.rast, p.geom)
>>>>
>>>> This returns a raster which has all pixels inside the polygon
>>>>
>>>> On Thu, Feb 5, 2015 at 4:05 PM, Jean Marchal <jean.d.marchal@gmail.com>
>>>> wrote:
>>>>
>>>>> Hi list,
>>>>>
>>>>> I am trying to return all the pixels in a raster that intersect (not
>>>>> just touch) an extent (say a rectangle). I tried ST_Clip and
>>>>> ST_Intersection(raster, geom) but they don't return all the pixels that
>>>>> intersect my extent polygon. Do I have to vectorize the raster first using
>>>>> ST_PixelAsPolygons or there is a better / more efficient way to proceed?
>>>>>
>>>>> Ultimately the goal is to fetch the resulting raster in R.
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Jean
>>>>>
>>>>> _______________________________________________
>>>>> postgis-users mailing list
>>>>> postgis-users@lists.osgeo.org
>>>>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> postgis-users mailing list
>>>> postgis-users@lists.osgeo.org
>>>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>>>
>>>
>>>
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users@lists.osgeo.org
>>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>>
>>
>>
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users@lists.osgeo.org
>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>
>
>
[Attachment #5 (text/html)]
<div dir="ltr"><div class="gmail_extra"><div><div>For the posterity, I modified very \
slightly the query by changing the returnband argument of ST_Intersection to return \
only one band.<br></div><div><br>WITH foo AS (</div><div><span \
style="white-space:pre-wrap"> </span>SELECT</div><div><span \
style="white-space:pre-wrap"> </span>ST_AsRaster(</div><span><div><span \
style="white-space:pre-wrap"> </span>ST_GeomFromText('POLYGON \
((-52.54178994517749 46.99199259385565, -52.54178994517749 46.996897959881, \
-52.53436080387823 46.996897959881, -52.53436080387823 46.99199259385565, \
-52.54178994517749 46.99199259385565))', 4269),</div></span><div><span \
style="white-space:pre-wrap"> </span>rast,</div><div><span \
style="white-space:pre-wrap"> </span>'8BUI',</div><div><span \
style="white-space:pre-wrap"> </span>touched := True</div><div><span \
style="white-space:pre-wrap"> </span>) AS rast</div><div><span \
style="white-space:pre-wrap"> </span>FROM elev</div><div><span \
style="white-space:pre-wrap"> </span>LIMIT \
1</div><div>)</div><div>SELECT</div><div><span \
style="white-space:pre-wrap"> </span>ST_Intersection(</div><div><span \
style="white-space:pre-wrap"> </span>elev.rast,</div><div><span \
style="white-space:pre-wrap"> </span>foo.rast,<br></div><div> \
'BAND1'<br></div><div><span style="white-space:pre-wrap"> </span>) AS \
rast</div><div>FROM elev</div><div>JOIN foo</div><div><span \
style="white-space:pre-wrap"> </span>ON ST_Intersects(elev.rast, \
foo.rast)<br><br></div><div>Thanks \
again,<br></div><div><br></div><div>Jean<br></div><div><br></div></div><div \
class="gmail_quote">2015-02-05 15:29 GMT-08:00 Jean Marchal <span dir="ltr"><<a \
href="mailto:jean.d.marchal@gmail.com" \
target="_blank">jean.d.marchal@gmail.com</a>></span>:<br><blockquote \
class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid \
rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>It worked! Thanks a \
lot!<br></div><div><br></div>Jean<br><div><div><div><br><div class="gmail_extra"><div \
class="gmail_quote">2015-02-05 15:18 GMT-08:00 Bborie Park <span dir="ltr"><<a \
href="mailto:dustymugs@gmail.com" \
target="_blank">dustymugs@gmail.com</a>></span>:<div><div><br><blockquote \
class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid \
rgb(204,204,204);padding-left:1ex"><div dir="ltr">Try something \
like:<br><div><br></div><div><div>WITH foo AS (</div><div><span \
style="white-space:pre-wrap"> </span>SELECT</div><div><span \
style="white-space:pre-wrap"> </span>ST_AsRaster(</div><span><div><span \
style="white-space:pre-wrap"> </span>ST_GeomFromText('POLYGON \
((-52.54178994517749 46.99199259385565, -52.54178994517749 46.996897959881, \
-52.53436080387823 46.996897959881, -52.53436080387823 46.99199259385565, \
-52.54178994517749 46.99199259385565))', 4269),</div></span><div><span \
style="white-space:pre-wrap"> </span>rast,</div><div><span \
style="white-space:pre-wrap"> </span>'8BUI',</div><div><span \
style="white-space:pre-wrap"> </span>touched := True</div><div><span \
style="white-space:pre-wrap"> </span>) AS rast</div><div><span \
style="white-space:pre-wrap"> </span>FROM elev</div><div><span \
style="white-space:pre-wrap"> </span>LIMIT \
1</div><div>)</div><div>SELECT</div><div><span \
style="white-space:pre-wrap"> </span>ST_Intersection(</div><div><span \
style="white-space:pre-wrap"> </span>elev.rast,</div><div><span \
style="white-space:pre-wrap"> </span>foo.rast</div><div><span \
style="white-space:pre-wrap"> </span>) AS rast</div><div>FROM elev</div><div>JOIN \
foo</div><div><span style="white-space:pre-wrap"> </span>ON ST_Intersects(elev.rast, \
foo.rast)</div></div><div><br></div><div>The idea is to explicitly convert the \
geometry into a raster using one of the elev rasters as a reference. The \
"touched := True" slightly changes the behavior of the \
rasterization.</div><div><br></div><div>From the docs:</div><div><br></div><div><span \
style="color:rgb(46,46,46);font-family:"Lucida \
Grande",Verdana,Geneva,Arial,Helvetica,sans-serif;font-size:13px">The optional \
</span><code style="color:rgb(46,46,46);font-size:13px">touched</code><span \
style="color:rgb(46,46,46);font-family:"Lucida \
Grande",Verdana,Geneva,Arial,Helvetica,sans-serif;font-size:13px"> parameter \
defaults to false and maps to the GDAL ALL_TOUCHED rasterization option, which \
determines if pixels touched by lines or polygons will be \
burned.</span><br></div></div><div><div><div class="gmail_extra"><br><div \
class="gmail_quote">On Thu, Feb 5, 2015 at 2:45 PM, Jean Marchal <span \
dir="ltr"><<a href="mailto:jean.d.marchal@gmail.com" \
target="_blank">jean.d.marchal@gmail.com</a>></span> wrote:<br><blockquote \
class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid \
rgb(204,204,204);padding-left:1ex"><div dir="ltr">Bborie, <br><br>I am running \
PostGIS 2.1.5. Here is the output of the \
postgis_full_version():<br><div><br>"POSTGIS="2.1.5 r13152" \
GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23 September 2009" \
GDAL="GDAL 1.9.0, released 2011/12/29" LIBXML="2.8.0" \
LIBJSON="UNKNOWN" RASTER"<br><br></div><div>My query looks like \
this:<br><br>SELECT ST_Intersection(ST_GeomFromText('POLYGON ((-52.54178994517749 \
46.99199259385565, -52.54178994517749 46.996897959881, -52.53436080387823 \
46.996897959881, -52.53436080387823 46.99199259385565, -52.54178994517749 \
46.99199259385565))', 4269), rast) as rast<br>FROM \
elev<br></div><div><br></div><div>David,<br><br></div><div>This query does not \
returns all the polygons that intersect my polygon. I think it returns only those \
where the centroid is inside the polygon or is it covered based? (like 50% of the \
pixel intersect with the polygon).<br></div><div><br></div><div>Thanks for your rapid \
answers!<br></div><div><div class="gmail_extra"><br></div><div \
class="gmail_extra">Jean<br></div><div class="gmail_extra"><br><div \
class="gmail_quote">2015-02-05 14:31 GMT-08:00 David Haynes <span dir="ltr"><<a \
href="mailto:haynesd2@gmail.com" \
target="_blank">haynesd2@gmail.com</a>></span>:<div><div><br><blockquote \
class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid \
rgb(204,204,204);padding-left:1ex"><div dir="ltr">select ST_Clip(r.rast,p.geom) as \
rast<div>from polygon p inner join raster r on ST_intersects(r.rast, \
p.geom)</div><div><br></div><div>This returns a raster which has all pixels inside \
the polygon</div></div><div class="gmail_extra"><br><div class="gmail_quote"><span>On \
Thu, Feb 5, 2015 at 4:05 PM, Jean Marchal <span dir="ltr"><<a \
href="mailto:jean.d.marchal@gmail.com" \
target="_blank">jean.d.marchal@gmail.com</a>></span> wrote:<br></span><blockquote \
class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid \
rgb(204,204,204);padding-left:1ex"><div><div><div dir="ltr"><div>Hi \
list,<br></div><div><br>I am trying to return all the pixels in a raster that \
intersect (not just touch) an extent (say a rectangle). I tried ST_Clip and \
ST_Intersection(raster, geom) but they don't return all the pixels that intersect \
my extent polygon. Do I have to vectorize the raster first using ST_PixelAsPolygons \
or there is a better / more efficient way to proceed? <br><br></div><div>Ultimately \
the goal is to fetch the resulting raster in \
R.<br><br></div><div>Thanks,<br></div><div><br>Jean<br></div></div> \
<br></div></div><span>_______________________________________________<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><br></span></blockquote></div><br></div>
<br>_______________________________________________<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><br></blockquote></div></div></div><br></div></div></div>
<br>_______________________________________________<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><br></blockquote></div><br></div>
</div></div><br>_______________________________________________<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><br></blockquote></div></div></div><br></div></div></div></div></div>
</blockquote></div><br></div></div>
_______________________________________________
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic