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
>>>
>>>
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users@postgis.refractions.net
>>> http://postgis.refractions.net/mailman/listinfo/postgis-users
--
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
bkp...@ucdavis.edu
_______________________________________________
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
