[postgis-users] ST_AsRaster and gdal_translate
Guillaume Drolet DRF
DROGU2 at intranet.mrn.gouv
Thu Jul 26 08:11:56 PDT 2012
Many thanks for your help!! You really pulled me out of my downward spiral of frustration ;)
Both suggestions from Pierre and Bborie work but the results are not quite the same although similar.
Burned pixels generally overlap between both solutions although in some cases pixels get burned with one method but not with the other and vice versa..
In this case, it doesn't really matter if values are only 1 or no data because I only need to get a raster of rivers and lakes (= 1)
versus no water (= no data).
Again, I couldn't thank you enough!
Gui
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 <http://postgis.refractions.net/mailman/listinfo/postgis-users> [mailto:postgis-users-
/>>/ bounces at postgis.refractions.net <http://postgis.refractions.net/mailman/listinfo/postgis-users>] On Behalf Of Bborie Park
/>>/ Sent: Wednesday, July 25, 2012 4:23 PM
/>>/ To:postgis-users at postgis.refractions.net <http://postgis.refractions.net/mailman/listinfo/postgis-users>
/>>/ 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>
/>>>/ http://postgis.refractions.net/mailman/listinfo/postgis-users
/>>>/
/>>/
/>>/ --
/>>/ Bborie Park
/>>/ Programmer
/>>/ Center for Vectorborne Diseases
/>>/ UC Davis
/>>/ 530-752-8380
/>>/ bkpark at ucdavis.edu <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>
/>>/ 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>
/>/ http://postgis.refractions.net/mailman/listinfo/postgis-users
/>/
/
--
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
bkpark at ucdavis.edu <http://postgis.refractions.net/mailman/listinfo/postgis-users>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20120726/0e9dfc7d/attachment.html>
More information about the postgis-users
mailing list