<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>