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 <bkp...@ucdavis.edu> >>> Subject: Re: [postgis-users] Problem using ST_AsRaster >>> To: postgis-users@postgis.refractions.net >>> Message-ID: <4fe8a268.6050...@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. >>>> >
_______________________________________________ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users
