[postgis-users] ST_AsRaster and gdal_translate

[Bborie Park]

By the looks of it, there is no burning of values from the polygon
layer, as per the call to ST_AsRaster().

CREATE TABLE rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m AS
    SELECT ST_AsRaster(hyd.geom, 1000.0, 1000.0) AS rast
    FROM hydrology.bdtq_20k_hydro_so_gfsm_aea_subset AS hyd;

That would just burn 1 or 0 (NODATA).

-bborie

On 07/25/2012 01:40 PM, Pierre Racine wrote:
> If he does that he will not get the right value from his polygon layer...
> 
> The right way to create an exportable raster from a polygon table it to:
> 
> 1) Rasterize and align properly your geometries (adding the gridx and gridy parameters to ST_AsRaster and the pixeltype and value to burn the right value):
> 
> CREATE TABLE rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m AS
> SELECT ST_AsRaster(hyd.geom, 1000.0, 1000.0, 0.0, 0.0, '32BUI'::text, nameOfTheColumnContainingTheValueToBurn ) AS rast
> FROM hydrology.bdtq_20k_hydro_so_gfsm_aea_subset AS hyd;
> 
> 2) Union the table.
> 
> If your raster is very small you can do:
> 
> CREATE TABLE rast2export AS
> SELECT ST_Union(rast) rast
> FROM rasters.bdtq_20k_hydro_so_gfsm_aea_subset
> 
> If your raster is big, ST_Union is too slow. You have to use the method I describe in my blog 
> 
> http://geospatialelucubrations.blogspot.ca/2012/07/a-slow-yet-1000x-faster-alternative-to.html
> 
> 3) Then you can GDAL_translate using mode=1
> 
> Sorry it is so complicated...
> 
> In the future there should just be a ST_UnionToRaster() function, planned in the specifications, that rasterize all the geometries to a common raster and then export. So just two minimal steps.
> 
> Pierre
> 
> 
> 
>> -----Original Message-----
>> From: postgis-users-bounces at postgis.refractions.net [mailto:postgis-users-
>> bounces at postgis.refractions.net] On Behalf Of Bborie Park
>> Sent: Wednesday, July 25, 2012 4:23 PM
>> To: postgis-users at postgis.refractions.net
>> Subject: Re: [postgis-users] ST_AsRaster and gdal_translate
>>
>> For exporting, instead of using the rasters in
>> "bdtq_20k_hydro_so_gfsm_aea_subset_1000m", you may want to try
>> something
>> like...
>>
>> CREATE TABLE foo AS
>> SELECT 1 AS rid, ST_AsRaster((
>> 	SELECT
>> 		ST_Collect(geom)
>> 	FROM hydrology.bdtq_20k_hydro_so_gfsm_aea_subset
>> ), 1000.0, 1000.0
>> ) AS rast
>>
>> The above creates one raster using the collection of geometries, which
>> can then be exported with gdal_translate.
>>
>> -bborie
>>
>> On 07/25/2012 12:57 PM, Guillaume Drolet DRF wrote:
>>> Hi,
>>>
>>> I have been struggling with this since yesterday and now I think it is
>>> time for me to request some help from the community :
>>>
>>> I am trying to convert a geometry table (polygons) to a raster that will
>>> have 1000 m-pixels, and then export that raster to a TIF file using
>>> gdal_translate.
>>> I've tried loads of different queries but nothing seems to work... Here
>>> are examples of some the queries I tried:
>>>
>>> CREATE TABLE rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m AS
>>>     SELECT ST_AsRaster(hyd.geom, 1000.0, 1000.0) AS rast
>>>     FROM hydrology.bdtq_20k_hydro_so_gfsm_aea_subset AS hyd;
>>>
>>> ALTER TABLE rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m ADD
>> COLUMN
>>> rid serial PRIMARY KEY;
>>>
>>> CREATE INDEX
>>>     ON rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m
>>>     USING gist(st_convexhull(rast));
>>>
>>> SELECT AddRasterConstraints('rasters'::name,
>>>         'bdtq_20k_hydro_so_gfsm_aea_subset_1000m'::name,
>>>         'rast'::name);
>>>
>>> Then at the command line:
>>>
>>>> gdal_translate "PG:host=localhost port=5432 dbname=testspdb
>>> user=postgres password=*** schema=rasters
>>> table=bdtq_20k_hydro_so_gfsm_aea_subset_1000m mode=2"
>>> E:\temp\hydro1000m.tif
>>>
>>> I get this error:
>>>
>>> ERROR 1: Error, the table
>>> rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m contains tiles with
>>> different size, and irregular blocking is not supported yet
>>> GDALOpen failed - 1
>>> Error, the table rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m
>>> contains tiles
>>> with different size, and irregular blocking is not supported yet
>>>
>>> My problem seems to be related to the irregular blocking of my raster. I
>>> though I could create a raster with regular block and pixel size but I
>>> didn't figure how to do this yet. I tried to achieve this by creating an
>>> empty raster:
>>>
>>> CREATE TABLE rasters.test_rasterize (
>>>     rid serial PRIMARY KEY,
>>>     rast raster
>>> );
>>>
>>> INSERT INTO rasters.test_rasterize VALUES (
>>>     1,
>>>     ST_MakeEmptyRaster(45, 62, 1828552.00, 1422569.00, 1000.0, 1000.0,
>>> 0.0, 0.0, 3175)
>>> );
>>> -- 45 and 62 are the approximate numbers of 1000 m pixels in x and y
>>> that fit inside the extent of my geometry layer
>>> UPDATE rasters.test_rasterize
>>>     SET rast = ST_AsRaster(geom, 1000.0, 1000.0)
>>>         FROM hydrology.bdtq_20k_hydro_so_gfsm_aea_subset;
>>>
>>>
>>> When I call gdal_translate with rasters.test_rasterize, I get this:
>>>
>>> Input file size is 2, 4
>>> 0...10...20...30...40...50...60...70...80...90...100 - done.
>>>
>>> I'm getting desperate and will greatly appreciate any hints!
>>>
>>> If that can help, I'm using PostGIS 2.0.0 r9605 and GDAL 1.9.0
>>>
>>> Many thanks,
>>>
>>> Guillaume
>>>
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users at postgis.refractions.net
>>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>>
>>
>> --
>> Bborie Park
>> Programmer
>> Center for Vectorborne Diseases
>> UC Davis
>> 530-752-8380
>> bkpark at ucdavis.edu
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at postgis.refractions.net
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
> 

-- 
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
bkpark at ucdavis.edu