[postgis-users] ST_AsRaster and gdal_translate

[Pierre Racine]

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