[prev in list] [next in list] [prev in thread] [next in thread]
List: postgis-users
Subject: Re: [postgis-users] Problem using ST_AsRaster: SOLVED
From: Mark Wynter <mark () dimensionaledge ! com>
Date: 2012-06-28 1:21:10
Message-ID: DE59CF1F-46A8-43BC-B1F7-D8DCCAADAA6D () dimensionaledge ! com
[Download RAW message or body]
[Attachment #2 (multipart/alternative)]
Thanks Bborie. I stuck with variant 1 but added an additional step of using \
ST_MapAlgebraExpr to burn the values to the udpated reference tile.
CREATE TABLE dummy_rast_updated AS
SELECT rt.rid, ST_MapAlgebraExpr(rt.rast, vtr.rast, '[rast2.val]', '8BUI', 'FIRST', \
'0', '0', '0') as rast FROM dummy_rast as rt, viewshed_rast as vtr;
The above expression gives the result I need. I can now string all this together \
as a 1-step process using pl/pgsql. :-)
>
> You could try variants 9 or 10 (the last two) described in the docs...
>
> http://postgis.refractions.net/documentation/manual-svn/RT_ST_AsRaster.html
>
> You'll probably need to call ST_Rescale afterwards to change the scale
> of the resulting raster to match that of your reference raster.
>
> By default, ST_AsRaster returns a raster having the extent of the input
> geometry's bounding box.
>
> -bborie
>
> On 06/25/2012 06:48 PM, Mark Wynter wrote:
> > Hi Bborie
> >
> > As part of a pl/pgsql update process, my requirement is to create a "new" raster \
> > tile with same alignment (AND upperleftx, upperlefty, block dimensions) as the \
> > reference raster tile. For the new tile, I then want to set the pixel values = \
> > 100 for those pixels which intersect the "updated" vector polygon. I'm not \
> > interested in parts of the vector polygon beyond the extent of the reference \
> > tile.
> > The reference raster is:
> > CREATE TABLE dummy_rast (rid integer, rast raster) WITH (OIDS=FALSE);
> > INSERT INTO dummy_rast VALUES(1, ST_MakeEmptyRaster(30, 30, 576000, -3780375, 50, \
> > -30, 0, 0, 3577)); UPDATE dummy_rast SET rast = ST_AddBand(rast, 1, '8BUI', \
> > NULL); SELECT AddRasterConstraints('dummy_rast'::name, 'rast'::name);
> >
> > As an approach, my initial thoughts were to rasterise the vector layer using \
> > ST_AsRaster.
> > From scratch, the vector layer comprises two overlapping polygons which protrude \
> > beyond the extent of the reference tile... CREATE TABLE viewshed_vectors (gid \
> > integer, geometry geometry(polygon, 3577)) WITH (OIDS=FALSE); INSERT INTO \
> > viewshed_vectors VALUES(1, \
> > ST_Polygon(ST_LineFromWKB(ST_AsBinary(ST_GeomFromText('LINESTRING(576500 \
> > -3780500, 577000 -3780500, 577000 -3781000, 576500 -3781000, 576500 \
> > -3780500)')),3577),3577)); INSERT INTO viewshed_vectors VALUES(2, \
> > ST_Polygon(ST_LineFromWKB(ST_AsBinary(ST_GeomFromText('LINESTRING(576750 \
> > -3780250, 577250 -3780250, 577250 -3780750, 576750 -3780750, 576750 \
> > -3780250)')),3577),3577));
> > To rasterise the viewshed_vectors...
> > CREATE TABLE viewshed_rast AS
> > WITH vt as (SELECT ST_Union(geometry) as geometry FROM viewshed_vectors)
> > SELECT rt.rid, ST_AsRaster(vt.geometry, rt.rast, '8BUI', 100, NULL) as rast FROM \
> > vt, dummy_rast as rt; SELECT AddRasterConstraints('viewshed_rast'::name, \
> > 'rast'::name);
> > SELECT ST_SameAlignment(rt.rast, vt.rast) as sm FROM dummy_rast as rt, \
> > viewshed_rast as vt; sm
> > ----
> > t
> >
> > The result when mapped in QGIS is as per this updated screenshot.
> >
> > Is there another step I must now do using ST_MapAlgebra to burn the \
> > (intersecting) viewshed_rast values to the underlying dummy_rast tile? What \
> > might that expression look like?
> > Should I be approaching this a different way?
> >
> > I appreciate your thoughts and help.
> >
> > Best regards
> >
> >
> >
> >
> > >
> > > Message: 18
> > > Date: Mon, 25 Jun 2012 10:39:52 -0700
> > > From: Bborie Park <bkpark@ucdavis.edu>
> > > Subject: Re: [postgis-users] Problem using ST_AsRaster
> > > To: postgis-users@postgis.refractions.net
> > > Message-ID: <4FE8A268.6050802@ucdavis.edu>
> > > Content-Type: text/plain; charset=ISO-8859-1
> > >
> > > Hi Mark,
> > >
> > > Can you elaborate on what you mean by "resultant raster does not map to
> > > the reference raster"?
> > >
> > > The output from ST_AsRaster should result in a raster with the same
> > > SRID, scale and skew as the reference raster. The output raster should
> > > also be aligned with the reference raster as tested by ST_SameAlignment.
> > >
> > > -bborie
> > >
> > > On 06/23/2012 11:04 PM, Mark Wynter wrote:
> > > > I can rasterise a vector layer, but I'm having trouble mapping it to a \
> > > > reference raster.
> > > > The reference raster, called dummy_rast is a 1x1 raster tile with a height \
> > > > and width of 500pixels, each of 250m in size. I created using a pl/pgsql \
> > > > function: SELECT make_tiled_raster('public', 'dummy_rast', 576000, -3780000, \
> > > > 1, 1, 500, 500, 250, -250); The result is
> > > >
> > > > srid | scale_x | scale_y | blocksize_x | blocksize_y | num_bands | \
> > > > pixel_types | nodata_values
> > > > ------+---------+---------+-------------+-------------+-----------+-------------+---------------
> > > > 3577 | 250 | -250 | 500 | 500 | 1 | {8BUI} \
> > > > | {NULL}
> > > >
> > > > I now wish to burn a vector layer onto this raster:
> > > >
> > > > CREATE TABLE viewshed_rast AS
> > > > WITH vt as (SELECT ST_Union(geometry) as geometry FROM viewshed_vectors)
> > > > SELECT rt.rid, ST_AsRaster(vt.geometry, rt.rast, '8BUI', 120, 100) as rast \
> > > > FROM dummy_rast as rt, vt;
> > > > The result is
> > > > srid | scale_x | scale_y | blocksize_x | blocksize_y | num_bands | \
> > > > pixel_types | nodata_values
> > > > ------+---------+---------+-------------+-------------+-----------+-------------+---------------
> > > > 3577 | 250 | -250 | 67 | 38 | 1 | {8BUI} \
> > > > | {100} (1 row)
> > > >
> > > > I do not understand why the resultant raster does not map to the reference \
> > > > raster? Refer screenshot attached showing the resultant layers in QGIS. \
> > > > The upperleftx and upperlefty, and the block size of the resultant raster are \
> > > > defined by the extent of the vector layer and not the reference raster.
> > > > Is there something obvious I'm doing wrong? Thanks.
> > > >
>
[Attachment #5 (unknown)]
<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; \
-webkit-line-break: after-white-space; "><div><div>Thanks Bborie. I stuck \
with variant 1 but added an additional step of using ST_MapAlgebraExpr to burn the \
values to the udpated reference tile. </div><div><br></div><div>CREATE TABLE \
dummy_rast_updated AS</div><div>SELECT rt.rid, ST_MapAlgebraExpr(rt.rast, vtr.rast, \
'[rast2.val]', '8BUI', 'FIRST', '0', '0', '0') as rast FROM dummy_rast as rt, \
viewshed_rast as vtr;</div><div><br></div><div>The above expression gives the result \
I need. I can now string all this together as a 1-step process using \
pl/pgsql. :-)</div><div><br></div><div><br></div></div><div><blockquote \
type="cite"><div><br>You could try variants 9 or 10 (the last two) described in the \
docs...<br><br><a href="http://postgis.refractions.net/documentation/manual-svn/RT_ST_ \
AsRaster.html">http://postgis.refractions.net/documentation/manual-svn/RT_ST_AsRaster.html</a><br><br>You'll \
probably need to call ST_Rescale afterwards to change the scale<br>of the resulting \
raster to match that of your reference raster.<br><br>By default, ST_AsRaster returns \
a raster having the extent of the input<br>geometry's bounding \
box.<br><br>-bborie<br><br>On 06/25/2012 06:48 PM, Mark Wynter wrote:<br><blockquote \
type="cite">Hi Bborie<br></blockquote><blockquote \
type="cite"><br></blockquote><blockquote type="cite">As part of a pl/pgsql update \
process, my requirement is to create a "new" raster tile with same alignment (AND \
upperleftx, upperlefty, block dimensions) as the reference raster tile. For the \
new tile, I then want to set the pixel values = 100 for those pixels \
which intersect the "updated" vector polygon. I'm not interested in \
parts of the vector polygon beyond the extent of the reference \
tile.<br></blockquote><blockquote type="cite"><br></blockquote><blockquote \
type="cite">The reference raster is:<br></blockquote><blockquote type="cite">CREATE \
TABLE dummy_rast (rid integer, rast raster) WITH \
(OIDS=FALSE);<br></blockquote><blockquote type="cite">INSERT INTO dummy_rast \
VALUES(1, ST_MakeEmptyRaster(30, 30, 576000, -3780375, 50, -30, 0, 0, \
3577));<br></blockquote><blockquote type="cite">UPDATE dummy_rast SET rast = \
ST_AddBand(rast, 1, '8BUI', NULL);<br></blockquote><blockquote \
type="cite">SELECT AddRasterConstraints('dummy_rast'::name, \
'rast'::name);<br></blockquote><blockquote type="cite"><br></blockquote><blockquote \
type="cite">As an approach, my initial thoughts were to rasterise the vector layer \
using ST_AsRaster.<br></blockquote><blockquote \
type="cite"><br></blockquote><blockquote type="cite">From scratch, the vector layer \
comprises two overlapping polygons which protrude beyond the extent of the reference \
tile...<br></blockquote><blockquote type="cite">CREATE TABLE viewshed_vectors (gid \
integer, geometry geometry(polygon, 3577)) WITH \
(OIDS=FALSE);<br></blockquote><blockquote type="cite">INSERT INTO viewshed_vectors \
VALUES(1, ST_Polygon(ST_LineFromWKB(ST_AsBinary(ST_GeomFromText('LINESTRING(576500 \
-3780500, 577000 -3780500, 577000 -3781000, 576500 -3781000, 576500 \
-3780500)')),3577),3577));<br></blockquote><blockquote type="cite">INSERT INTO \
viewshed_vectors VALUES(2, \
ST_Polygon(ST_LineFromWKB(ST_AsBinary(ST_GeomFromText('LINESTRING(576750 -3780250, \
577250 -3780250, 577250 -3780750, 576750 -3780750, 576750 \
-3780250)')),3577),3577));<br></blockquote><blockquote \
type="cite"><br></blockquote><blockquote type="cite">To rasterise the \
viewshed_vectors...<br></blockquote><blockquote type="cite">CREATE TABLE \
viewshed_rast AS<br></blockquote><blockquote type="cite">WITH vt as (SELECT \
ST_Union(geometry) as geometry FROM viewshed_vectors)<br></blockquote><blockquote \
type="cite">SELECT rt.rid, ST_AsRaster(vt.geometry, rt.rast, '8BUI', 100, NULL) as \
rast FROM vt, dummy_rast as rt;<br></blockquote><blockquote type="cite">SELECT \
AddRasterConstraints('viewshed_rast'::name, \
'rast'::name);<br></blockquote><blockquote type="cite"><br></blockquote><blockquote \
type="cite">SELECT ST_SameAlignment(rt.rast, vt.rast) as sm FROM dummy_rast as rt, \
viewshed_rast as vt;<br></blockquote><blockquote type="cite">sm \
<br></blockquote><blockquote type="cite">----<br></blockquote><blockquote \
type="cite"> t<br></blockquote><blockquote type="cite"><br></blockquote><blockquote \
type="cite">The result when mapped in QGIS is as per this updated screenshot. \
<br></blockquote><blockquote \
type="cite"><br></blockquote><blockquote type="cite">Is there another step I must now \
do using ST_MapAlgebra to burn the (intersecting) viewshed_rast values to the \
underlying dummy_rast tile? What might that expression look like? \
<br></blockquote><blockquote type="cite"><br></blockquote><blockquote \
type="cite">Should I be approaching this a different way? \
<br></blockquote><blockquote \
type="cite"><br></blockquote><blockquote type="cite">I appreciate your thoughts and \
help.<br></blockquote><blockquote type="cite"><br></blockquote><blockquote \
type="cite">Best regards<br></blockquote><blockquote \
type="cite"><br></blockquote><blockquote type="cite"><br></blockquote><blockquote \
type="cite"><br></blockquote><blockquote type="cite"><br></blockquote><blockquote \
type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite">Message: \
18<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">Date: \
Mon, 25 Jun 2012 10:39:52 -0700<br></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite">From: Bborie Park \
<bkpark@ucdavis.edu><br></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite">Subject: Re: [postgis-users] Problem using \
ST_AsRaster<br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite">To: postgis-users@postgis.refractions.net<br></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite">Message-ID: \
<4FE8A268.6050802@ucdavis.edu><br></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite">Content-Type: text/plain; \
charset=ISO-8859-1<br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite">Hi Mark,<br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite">Can you elaborate on what you mean by "resultant raster does not map \
to<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">the \
reference raster"?<br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite">The output from ST_AsRaster should result in a raster with the \
same<br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite">SRID, scale and skew as the reference raster. The output raster \
should<br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite">also be aligned with the reference raster as tested by \
ST_SameAlignment.<br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite">-bborie<br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite">On 06/23/2012 11:04 PM, Mark Wynter \
wrote:<br></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote type="cite">I can rasterise a vector layer, but I'm having \
trouble mapping it to a reference \
raster.<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote \
type="cite"><br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">The reference raster, \
called dummy_rast is a 1x1 raster tile with a height and width of 500pixels, each of \
250m in size. I created using a pl/pgsql \
function:<br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">SELECT \
make_tiled_raster('public', 'dummy_rast', 576000, -3780000, 1, 1, 500, 500, 250, \
-250);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote type="cite">The result is \
<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote \
type="cite"><br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">srid | scale_x | scale_y \
| blocksize_x | blocksize_y | num_bands | pixel_types | nodata_values \
<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote \
type="cite">------+---------+---------+-------------+-------------+-----------+-------------+---------------<br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">3577 | \
250 | -250 | \
500 | \
500 | \
1 | {8BUI} \
| \
{NULL}<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote \
type="cite"><br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote \
type="cite"><br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">I now wish to burn a \
vector layer onto this raster:<br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote \
type="cite"><br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">CREATE TABLE \
viewshed_rast AS<br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">WITH vt as (SELECT \
ST_Union(geometry) as geometry FROM \
viewshed_vectors)<br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">SELECT rt.rid, \
ST_AsRaster(vt.geometry, rt.rast, '8BUI', 120, 100) as rast FROM dummy_rast as rt, \
vt;<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote \
type="cite"><br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">The result is \
<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote type="cite">srid | scale_x | scale_y | blocksize_x | \
blocksize_y | num_bands | pixel_types | nodata_values \
<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote \
type="cite">------+---------+---------+-------------+-------------+-----------+-------------+---------------<br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">3577 | \
250 | -250 | \
67 | \
38 | \
1 | {8BUI} \
| \
{100}<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote type="cite">(1 \
row)<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote \
type="cite"><blockquote \
type="cite"><br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">I do not understand why \
the resultant raster does not map to the reference raster? Refer \
screenshot attached showing the resultant layers in QGIS. The upperleftx and \
upperlefty, and the block size of the resultant raster are defined by the extent of \
the vector layer and not the reference raster. \
<br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote \
type="cite"><br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote type="cite">Is there something \
obvious I'm doing wrong? \
Thanks.<br></blockquote></blockquote></blockquote><blockquote \
type="cite"><blockquote type="cite"><blockquote \
type="cite"><br></blockquote></blockquote></blockquote><font class="Apple-style-span" \
color="#540000"><br></font></div></blockquote></div><br></body></html>
_______________________________________________
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/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