<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"><head><meta http-equiv=Content-Type content="text/html; charset=iso-8859-1"><meta name=Generator content="Microsoft Word 15 (filtered medium)"><style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:#0563C1;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:#954F72;
text-decoration:underline;}
span.EmailStyle17
{mso-style-type:personal;
font-family:"Calibri",sans-serif;
color:windowtext;}
span.EmailStyle18
{mso-style-type:personal;
font-family:"Calibri",sans-serif;
color:#1F497D;}
span.EmailStyle20
{mso-style-type:personal-reply;
font-family:"Calibri",sans-serif;
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:8.5in 11.0in;
margin:70.85pt 85.05pt 70.85pt 85.05pt;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]--></head><body lang=EN-US link="#0563C1" vlink="#954F72"><div class=WordSection1><p class=MsoNormal><span style='color:#1F497D'>Just occurred to me that the ST_Intersection I gave you earlier probably does much the same as ST_Clip.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>That said, you should be able to achieve a less choppy result using ST_Resample and ST_Clip in unison.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Increase the 5 to as much as you need to get the less choppy effect you are looking for.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal>SELECT ST_Clip( <span style='color:#1F497D'>ST_Resample(rast, ST_Width(rast)*5, ST_Height(rast)*5 ) , c.</span> area_country<span style='color:#1F497D'>) <o:p></o:p></span></p><p class=MsoNormal style='text-indent:.5in'>FROM country AS c<o:p></o:p></p><p class=MsoNormal style='margin-left:.5in;text-indent:.5in'>INNER JOIN tb_densitysurface AS d <o:p></o:p></p><p class=MsoNormal style='margin-left:1.0in;text-indent:.5in'>ON ( ST_Intersects(d.rast, c.area_country) )<o:p></o:p></p><p class=MsoNormal>WHERE c.gid = 19<span lang=PT-BR>;</span><o:p></o:p></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><div style='border:none;border-left:solid blue 1.5pt;padding:0in 0in 0in 4.0pt'><div><div style='border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in'><p class=MsoNormal><b>From:</b> Regina Obe [mailto:lr@pcorp.us] <br><b>Sent:</b> Sunday, December 4, 2022 12:10 PM<br><b>To:</b> 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org><br><b>Cc:</b> 'liglio.pessoal@nexxa.com.br' <liglio.pessoal@nexxa.com.br><br><b>Subject:</b> RE: [postgis-users] RES: Clipping a raster and a polygon using ST_Intersection<o:p></o:p></p></div></div><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><span style='color:#1F497D'>Liglio,<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>b1 being choppier than b2, is because the rasterized polygon takes on the pixel size of your low res polygon.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>I’d try resampling your </span>tb_densitysurface to lower pixel size before you do the operation with ST_Resample.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><span style='color:#1F497D'><a href="https://postgis.net/docs/RT_ST_Resample.html">https://postgis.net/docs/RT_ST_Resample.html</a><o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>The width and height are measured in pixels. The idea is you want your tb_densitysurface to have more pixels than it does now so that the size of each pixel is smaller.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>If you do something like <o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>CREATE TABLE tb_densitysurface_x5 AS <o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>SELECT rid, ST_Resample(rast, ST_Width(rast)*5, ST_Height(rast)*5 ) AS rast<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>FROM tb_densitysurface;<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>That should give you the same spatial ref sys, but the pixels will be smaller in size.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>You can then use this to rerun the ST_Intersection in raster space. <o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>I think it will be a little slower though, because the more pixels you have the slower things get.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Just keep on increasing that number 5 until you get something that meets your choppiness criteria.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>As you get less choppy, operations will get slower, so you need to find a sweet spot there.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Regarding: </span>There is a way to export this image and associate a value for each square (polygon) to be colored accordingly ?<o:p></o:p></p><p class=MsoNormal>Was that to resolve the above issue, or you just want to know.<o:p></o:p></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Well there is ST_PixelAsPolygons which will convert each pixel to a square polygon. <a href="https://postgis.net/docs/RT_ST_PixelAsPolygons.html">https://postgis.net/docs/RT_ST_PixelAsPolygons.html</a><o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>There is also ST_DumpAsPolygons <a href="https://postgis.net/docs/RT_ST_DumpAsPolygons.html">https://postgis.net/docs/RT_ST_DumpAsPolygons.html</a> which tries to aggregate nearby pixels with the same value into a polygon.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Both return geomval rows similar to the ST_Intersection of raster, geom. As I recall I think those were completely rewritten in C so should still be faster than ST_Intersection(geometry,raster).<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><div style='border:none;border-left:solid blue 1.5pt;padding:0in 0in 0in 4.0pt'><div><div style='border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in'><p class=MsoNormal><b>From:</b> postgis-users [<a href="mailto:postgis-users-bounces@lists.osgeo.org">mailto:postgis-users-bounces@lists.osgeo.org</a>] <b>On Behalf Of </b><a href="mailto:liglio.pessoal@nexxa.com.br">liglio.pessoal@nexxa.com.br</a><br><b>Sent:</b> Wednesday, November 30, 2022 9:37 AM<br><b>To:</b> 'PostGIS Users Discussion' <<a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a>><br><b>Subject:</b> [postgis-users] RES: Clipping a raster and a polygon using ST_Intersection<o:p></o:p></p></div></div><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><span lang=PT-BR>Regina,<o:p></o:p></span></p><p class=MsoNormal><span lang=PT-BR><o:p> </o:p></span></p><p class=MsoNormal>Thanks for your replay. <o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Your solution is fast like using ST_Clip, but I still have the problem that my pixels are 25km square (low definition) and the image look like this: (image b1.jpg)<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>And I want to produce something like this: (image b2.jpg)<o:p></o:p></p><p class=MsoNormal><b><o:p> </o:p></b></p><p class=MsoNormal>I know that the second image cannot be reproduced as just raster with big pixels, the pixels cannot be cutted.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>There is a way to export this image and associate a value for each square (polygon) to be colored accordingly ?<o:p></o:p></p><div><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Regards,<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Liglio<o:p></o:p></p></div><p class=MsoNormal><o:p> </o:p></p><div><div style='border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in'><p class=MsoNormal><b><span lang=PT-BR>De:</span></b><span lang=PT-BR> postgis-users <<a href="mailto:postgis-users-bounces@lists.osgeo.org">postgis-users-bounces@lists.osgeo.org</a>> <b>Em nome de </b>Regina Obe<br><b>Enviada em:</b> quinta-feira, 24 de novembro de 2022 12:19<br><b>Para:</b> 'PostGIS Users Discussion' <<a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a>><br><b>Assunto:</b> Re: [postgis-users] Clipping a raster and a polygon using ST_Intersection<o:p></o:p></span></p></div></div><p class=MsoNormal><span lang=PT-BR><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>I’m not quite sure what you mean by edges aren’t smooth. Rasters are never really smooth. Vectors are smooth. Rasters may look smooth at the edges, but that’s more because pixels are square and the pixels in each tile line up with the pixels.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>If you want it to look a little smoother, maybe expand your envelop ever so slightly so it lines up with a pixel in a raster. So each edge pixel is guaranteed to be fully inside the raster and no edge partly outside the raster or do the operation in raster space which I will explain shortly.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>ST_Intersection yes is known to be very slow and as you noticed does not give you a raster back. <o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Is ST_Clip + ST_Intersects performance sufficient for your needs? <o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>ST_Clip really is ST_Intersection in raster space for raster/geometry. So you could do the intersection in pure raster space by converting your geometry to a raster (making sure to align your new raster with your reference raster). That should make the edges smooth.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Like so – not tested so might have some syntax errors and also I haven’t used this in a very long time so may be wrong about what it returns.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal> SELECT ST_Intersection(r.rast, r2.rast) AS new_rast<o:p></o:p></p><p class=MsoNormal> FROM tb_densitysurface AS r, country, <o:p></o:p></p><p class=MsoNormal>-- create a raster to represent the country that has the same alignment as our density raster (pixels use the same grid)<o:p></o:p></p><p class=MsoNormal> ST_AsRaster( country.area_country, r.rast, <var><span style='font-family:"Calibri",sans-serif'>'1BB', 0</span></var>) AS r2(rast)<o:p></o:p></p><p class=MsoNormal> WHERE country.gid = 19 AND ST_Intersects(r.rast, r2.rast) <o:p></o:p></p><p class=MsoNormal><span lang=PT-BR>) As foo<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>As far as tiling goes, yah that would improve the smoothness of the inner tiles but not most outer ones.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>You’d also need to ST_Union them all up to get back to a single tile which will add perhaps more overhead than just keeping as single tile.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>If you want to tile your raster, you can use this function - <a href="https://postgis.net/docs/manual-2.5/RT_ST_Tile.html">https://postgis.net/docs/manual-2.5/RT_ST_Tile.html</a><o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Hope that helps,<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Regina<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='color:#1F497D'><o:p> </o:p></span></p><div style='border:none;border-left:solid blue 1.5pt;padding:0in 0in 0in 4.0pt'><div><div style='border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in'><p class=MsoNormal><b>From:</b> postgis-users [<a href="mailto:postgis-users-bounces@lists.osgeo.org">mailto:postgis-users-bounces@lists.osgeo.org</a>] <b>On Behalf Of </b><a href="mailto:liglio.pessoal@nexxa.com.br">liglio.pessoal@nexxa.com.br</a><br><b>Sent:</b> Tuesday, November 22, 2022 2:00 PM<br><b>To:</b> <a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br><b>Subject:</b> [postgis-users] Clipping a raster and a polygon using ST_Intersection<o:p></o:p></p></div></div><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Hi,<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>I created a retangular raster representing number of lightnings in a square. So I have this matrix raster created using ST_MapAlgebraFct and ST_MakeEmptyRaster. <o:p></o:p></p><p class=MsoNormal>This raster is not a tiled raster (tb_densitysurface), and when I Clip this raster using a polygon, the query below, the performance is not good.<o:p></o:p></p><p class=MsoNormal>Maybe because the raster is not tiled. So how do I tile this raster ? And I want the result to be a raster (not gval) to be exported as a file, and then colored using GDAL.<o:p></o:p></p><p class=MsoNormal>It worked with a different query only using ST_Clip (this returns a rast column and is very fast), but the edges arenīt smooth as using ST_Insersection. I am using Postgis 2.5.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>SELECT<o:p></o:p></p><p class=MsoNormal> (gval).val<o:p></o:p></p><p class=MsoNormal> (gval).geom<o:p></o:p></p><p class=MsoNormal>FROM (<o:p></o:p></p><p class=MsoNormal> SELECT ST_Intersection(<o:p></o:p></p><p class=MsoNormal> ST_Clip(rast,ST_Envelope(area_country)),<o:p></o:p></p><p class=MsoNormal> 1,<o:p></o:p></p><p class=MsoNormal> area_country<o:p></o:p></p><p class=MsoNormal> ) As gval<o:p></o:p></p><p class=MsoNormal> FROM country, tb_densitysurface WHERE ST_Intersects(rast, area_country) and country.gid = 19<o:p></o:p></p><p class=MsoNormal><span lang=PT-BR>) As foo<o:p></o:p></span></p><p class=MsoNormal><span lang=PT-BR><o:p> </o:p></span></p><p class=MsoNormal><span lang=PT-BR><o:p> </o:p></span></p><p class=MsoNormal><span lang=PT-BR>Regards,<o:p></o:p></span></p><p class=MsoNormal><span lang=PT-BR><o:p> </o:p></span></p><p class=MsoNormal><span lang=PT-BR>Liglio<o:p></o:p></span></p><p class=MsoNormal><span lang=PT-BR><o:p> </o:p></span></p></div></div></div></div></body></html>