[prev in list] [next in list] [prev in thread] [next in thread]
List: postgis-users
Subject: Re: [postgis-users] RES: Clipping a raster and a polygon using ST_Intersection
From: "Regina Obe" <lr () pcorp ! us>
Date: 2022-12-06 19:56:05
Message-ID: 000d01d909ad$c1f28930$45d79b90$ () pcorp ! us
[Download RAW message or body]
This is a multipart message in MIME format.
[Attachment #2 (multipart/alternative)]
This is a multipart message in MIME format.
ST_Resize might be a better option
https://postgis.net/docs/RT_ST_Resize.html
ST_Resize takes in raster width/height similar to ST_Resample.
From: Regina Obe [mailto:lr@pcorp.us]
Sent: Tuesday, December 6, 2022 2:53 PM
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Cc: 'liglio.pessoal@nexxa.com.br' <liglio.pessoal@nexxa.com.br>
Subject: RE: [postgis-users] RES: Clipping a raster and a polygon using
ST_Intersection
> Each pixel is associated with a value in the original raster, and when you
do resample, what value is associated to the pixels (in the docs “New pixel
values are computed using the NearestNeighbor”).
So in your example x5, only one pixel has the original pixel value and the
others 24 pixels has zero or nodata ? If I do colored scale of values it not
seems right.
Hmm that is not what I was expecting. Then again I always use the function
to reduce resolution, and not for this particular purpose. I had assumed it
would apply the value to the next neighboring and keep on applying that way.
But sounds like I was wrong about that.
I think I was confusing that with ST_Rescale. Try using ST_Rescale instead
and see if that gives expected results.
Note ST_Rescale takes as input pixel sizes. In your case you want to reduce
the size of each pixel. The below reduces the pixel size by 20%
https://postgis.net/docs/RT_ST_Rescale.html
SELECT ST_Clip( ST_Rescale(rast, ST_PixelWidth(rast)*0.2,
ST_PixelHeight(rast)*0.2 ) , c. area_country)
FROM country AS c
INNER JOIN tb_densitysurface AS d
ON ( ST_Intersects(d.rast, c.area_country) )
WHERE c.gid = 19;
Another question: I do a retangular matrix larger than the country (creating
an empty raster and dumping each pixel as polygon, and counting lightining
inside each sqare). So I have a lot of squares to compute, and 20 years of
lightining (about 2 billions).
It will be more efficient/faster if my empty raster has the shape of my
country, so It will have less squares to compute. It is possible ?
Yes it would be more efficient both in space and time for computing.
Just clip your raster by the country boundaries kind of like you were doing.
So something like
CREATE TABLE tb_densitysurface _country_clip AS
SELECT s.rid, ST_Clip(s.rast, geom) AS rast
FROM tb_densitysurface AS s CROSS JOIN (SELECT ST_Union(area_country) AS
geom FROM country_boundaries) AS c;
I recall you said you had just one raster, if you have multiple, then you
should change the above an INNER JOIN
ST_Intersects(s.rast, c.geom). If you only care about one country, you can
clip to that one country.
The only reason I am unioning is so you don’t end up with duplicate rasters
at country boundaries.
Hope that helps,
Regina
From: liglio.pessoal@nexxa.com.br <mailto:liglio.pessoal@nexxa.com.br>
[mailto:liglio.pessoal@nexxa.com.br]
Sent: Tuesday, December 6, 2022 2:20 PM
To: 'Regina Obe' <lr@pcorp.us <mailto:lr@pcorp.us> >
Subject: RES: [postgis-users] RES: Clipping a raster and a polygon using
ST_Intersection
Regina,
Thank you very much for your reply, and by the way congratulations for your
book “PostGis in action” third edition, I bought and it is very good.
Just occurred to me that the ST_Intersection I gave you earlier probably
does much the same as ST_Clip
R: Yes is pretty much the same. The result image is the same.
Each pixel is associated with a value in the original raster, and when you
do resample, what value is associated to the pixels (in the docs “New pixel
values are computed using the NearestNeighbor”).
So in your example x5, only one pixel has the original pixel value and the
others 24 pixels has zero or nodata ? If I do colored scale of values it not
seems right.
Another question: I do a retangular matrix larger than the country (creating
an empty raster and dumping each pixel as polygon, and counting lightining
inside each sqare). So I have a lot of squares to compute, and 20 years of
lightining (about 2 billions).
It will be more efficient/faster if my empty raster has the shape of my
country, so It will have less squares to compute. It is possible ?
Regards,
Liglio
De: Regina Obe <lr@pcorp.us <mailto:lr@pcorp.us> >
Enviada em: segunda-feira, 5 de dezembro de 2022 13:25
Para: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org
<mailto:postgis-users@lists.osgeo.org> >
Cc: liglio.pessoal@nexxa.com.br <mailto:liglio.pessoal@nexxa.com.br>
Assunto: RE: [postgis-users] RES: Clipping a raster and a polygon using
ST_Intersection
Just occurred to me that the ST_Intersection I gave you earlier probably
does much the same as ST_Clip.
That said, you should be able to achieve a less choppy result using
ST_Resample and ST_Clip in unison.
Increase the 5 to as much as you need to get the less choppy effect you are
looking for.
SELECT ST_Clip( ST_Resample(rast, ST_Width(rast)*5, ST_Height(rast)*5 ) , c.
area_country)
FROM country AS c
INNER JOIN tb_densitysurface AS d
ON ( ST_Intersects(d.rast, c.area_country) )
WHERE c.gid = 19;
From: Regina Obe [mailto:lr@pcorp.us]
Sent: Sunday, December 4, 2022 12:10 PM
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org
<mailto:postgis-users@lists.osgeo.org> >
Cc: 'liglio.pessoal@nexxa.com.br' <liglio.pessoal@nexxa.com.br
<mailto:liglio.pessoal@nexxa.com.br> >
Subject: RE: [postgis-users] RES: Clipping a raster and a polygon using
ST_Intersection
Liglio,
b1 being choppier than b2, is because the rasterized polygon takes on the
pixel size of your low res polygon.
I’d try resampling your tb_densitysurface to lower pixel size before you do
the operation with ST_Resample.
https://postgis.net/docs/RT_ST_Resample.html
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.
If you do something like
CREATE TABLE tb_densitysurface_x5 AS
SELECT rid, ST_Resample(rast, ST_Width(rast)*5, ST_Height(rast)*5 ) AS rast
FROM tb_densitysurface;
That should give you the same spatial ref sys, but the pixels will be
smaller in size.
You can then use this to rerun the ST_Intersection in raster space.
I think it will be a little slower though, because the more pixels you have
the slower things get.
Just keep on increasing that number 5 until you get something that meets
your choppiness criteria.
As you get less choppy, operations will get slower, so you need to find a
sweet spot there.
Regarding: There is a way to export this image and associate a value for
each square (polygon) to be colored accordingly ?
Was that to resolve the above issue, or you just want to know.
Well there is ST_PixelAsPolygons which will convert each pixel to a square
polygon. https://postgis.net/docs/RT_ST_PixelAsPolygons.html
There is also ST_DumpAsPolygons
https://postgis.net/docs/RT_ST_DumpAsPolygons.html which tries to aggregate
nearby pixels with the same value into a polygon.
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).
From: postgis-users [mailto:postgis-users-bounces@lists.osgeo.org] On Behalf
Of liglio.pessoal@nexxa.com.br <mailto:liglio.pessoal@nexxa.com.br>
Sent: Wednesday, November 30, 2022 9:37 AM
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org
<mailto:postgis-users@lists.osgeo.org> >
Subject: [postgis-users] RES: Clipping a raster and a polygon using
ST_Intersection
Regina,
Thanks for your replay.
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)
And I want to produce something like this: (image b2.jpg)
I know that the second image cannot be reproduced as just raster with big
pixels, the pixels cannot be cutted.
There is a way to export this image and associate a value for each square
(polygon) to be colored accordingly ?
Regards,
Liglio
De: postgis-users <postgis-users-bounces@lists.osgeo.org
<mailto:postgis-users-bounces@lists.osgeo.org> > Em nome de Regina Obe
Enviada em: quinta-feira, 24 de novembro de 2022 12:19
Para: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org
<mailto:postgis-users@lists.osgeo.org> >
Assunto: Re: [postgis-users] Clipping a raster and a polygon using
ST_Intersection
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.
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.
ST_Intersection yes is known to be very slow and as you noticed does not
give you a raster back.
Is ST_Clip + ST_Intersects performance sufficient for your needs?
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.
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.
SELECT ST_Intersection(r.rast, r2.rast) AS new_rast
FROM tb_densitysurface AS r, country,
-- create a raster to represent the country that has the same alignment as
our density raster (pixels use the same grid)
ST_AsRaster( country.area_country, r.rast, '1BB', 0) AS
r2(rast)
WHERE country.gid = 19 AND ST_Intersects(r.rast, r2.rast)
) As foo
As far as tiling goes, yah that would improve the smoothness of the inner
tiles but not most outer ones.
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.
If you want to tile your raster, you can use this function -
https://postgis.net/docs/manual-2.5/RT_ST_Tile.html
Hope that helps,
Regina
From: postgis-users [mailto:postgis-users-bounces@lists.osgeo.org] On Behalf
Of liglio.pessoal@nexxa.com.br <mailto:liglio.pessoal@nexxa.com.br>
Sent: Tuesday, November 22, 2022 2:00 PM
To: postgis-users@lists.osgeo.org <mailto:postgis-users@lists.osgeo.org>
Subject: [postgis-users] Clipping a raster and a polygon using
ST_Intersection
Hi,
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.
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.
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.
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.
SELECT
(gval).val
(gval).geom
FROM (
SELECT ST_Intersection(
ST_Clip(rast,ST_Envelope(area_country)),
1,
area_country
) As gval
FROM country, tb_densitysurface WHERE ST_Intersects(rast, area_country)
and country.gid = 19
) As foo
Regards,
Liglio
[Attachment #5 (text/html)]
<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: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'>ST_Resize might be a better option<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'><a \
href="https://postgis.net/docs/RT_ST_Resize.html">https://postgis.net/docs/RT_ST_Resize.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'>ST_Resize takes in raster width/height \
similar to ST_Resample.<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> \
Tuesday, December 6, 2022 2:53 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'>> </span>Each pixel is associated with a value in the \
original raster, and when you do resample, what value is associated to the pixels (in \
the docs “New pixel values are computed using the NearestNeighbor”). \
<o:p></o:p></p><p class=MsoNormal>So in your example x5, only one pixel has the \
original pixel value and the others 24 pixels has zero or nodata ? If I do colored \
scale of values it not seems right.<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'>Hmm that is not what I was expecting. Then again I always \
use the function to reduce resolution, and not for this particular purpose. I \
had assumed it would apply the value to the next neighboring and keep on applying \
that way. But sounds like I was wrong about that.<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'>I think I was confusing that with \
ST_Rescale. Try using ST_Rescale instead and see if that gives expected \
results.<o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Note \
ST_Rescale takes as input pixel sizes. In your case you want to reduce the size \
of each pixel. The below reduces the pixel size by 20%<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'><a \
href="https://postgis.net/docs/RT_ST_Rescale.html">https://postgis.net/docs/RT_ST_Rescale.html</a><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_Rescale(rast, \
ST_PixelWidth(rast)*0.2, ST_PixelHeight(rast)*0.2 ) , 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>Another \
question: I do a retangular matrix larger than the country (creating an empty raster \
and dumping each pixel as polygon, and counting lightining inside each sqare). So I \
have a lot of squares to compute, and 20 years of lightining (about 2 billions). \
<o:p></o:p></p><p class=MsoNormal>It will be more efficient/faster if my empty raster \
has the shape of my country, so It will have less squares to compute. It is possible \
?<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'>Yes it would be more efficient both in space and time for \
computing. <o:p></o:p></span></p><p class=MsoNormal><span style='color:#1F497D'>Just \
clip your raster by the country boundaries kind of like you were \
doing.<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'>So 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 </span>tb_densitysurface _country_clip AS \
<o:p></o:p></p><p class=MsoNormal>SELECT s.rid, ST_Clip(s.rast, geom) AS \
rast<o:p></o:p></p><p class=MsoNormal>FROM tb_densitysurface AS s CROSS JOIN \
(SELECT ST_Union(area_country) AS geom FROM country_boundaries) AS \
c;<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>I recall \
you said you had just one raster, if you have multiple, then you should change the \
above an INNER JOIN <o:p></o:p></p><p class=MsoNormal>ST_Intersects(s.rast, \
c.geom). If you only care about one country, you can clip to that one \
country.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>The \
only reason I am unioning is so you don’t end up with duplicate rasters at \
country boundaries.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p \
class=MsoNormal>Hope that helps,<o:p></o:p></p><p class=MsoNormal>Regina <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> <a \
href="mailto:liglio.pessoal@nexxa.com.br">liglio.pessoal@nexxa.com.br</a> [<a \
href="mailto:liglio.pessoal@nexxa.com.br">mailto:liglio.pessoal@nexxa.com.br</a>] \
<br><b>Sent:</b> Tuesday, December 6, 2022 2:20 PM<br><b>To:</b> 'Regina Obe' <<a \
href="mailto:lr@pcorp.us">lr@pcorp.us</a>><br><b>Subject:</b> RES: [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>Thank \
you very much for your reply, and by the way congratulations for your book \
“PostGis in action” third edition, I bought and it is very \
good.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><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</span><o:p></o:p></p><p class=MsoNormal>R: Yes \
is pretty much the same. The result image is the same.<o:p></o:p></p><p \
class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Each pixel is associated with \
a value in the original raster, and when you do resample, what value is associated to \
the pixels (in the docs “New pixel values are computed using the \
NearestNeighbor”). <o:p></o:p></p><p class=MsoNormal>So in your example x5, \
only one pixel has the original pixel value and the others 24 pixels has zero or \
nodata ? If I do colored scale of values it not seems right.<o:p></o:p></p><p \
class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Another question: I do a \
retangular matrix larger than the country (creating an empty raster and dumping each \
pixel as polygon, and counting lightining inside each sqare). So I have a lot of \
squares to compute, and 20 years of lightining (about 2 billions). <o:p></o:p></p><p \
class=MsoNormal>It will be more efficient/faster if my empty raster has the shape of \
my country, so It will have less squares to compute. It is possible \
?<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><div><p \
class=MsoNormal>Regards,<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p \
class=MsoNormal><span lang=PT-BR>Liglio<o:p></o:p></span></p></div><p \
class=MsoNormal><span lang=PT-BR><o:p> </o:p></span></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> Regina Obe <<a \
href="mailto:lr@pcorp.us">lr@pcorp.us</a>> <br><b>Enviada em:</b> segunda-feira, 5 \
de dezembro de 2022 13:25<br><b>Para:</b> 'PostGIS Users Discussion' <<a \
href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a>><br><b>Cc:</b> \
<a href="mailto:liglio.pessoal@nexxa.com.br">liglio.pessoal@nexxa.com.br</a><br><b>Assunto:</b> \
RE: [postgis-users] RES: 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'>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 [<a \
href="mailto:lr@pcorp.us">mailto:lr@pcorp.us</a>] <br><b>Sent:</b> Sunday, December \
4, 2022 12:10 PM<br><b>To:</b> 'PostGIS Users Discussion' <<a \
href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a>><br><b>Cc:</b> \
'liglio.pessoal@nexxa.com.br' <<a \
href="mailto:liglio.pessoal@nexxa.com.br">liglio.pessoal@nexxa.com.br</a>><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 \
_______________________________________________
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/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