[prev in list] [next in list] [prev in thread] [next in thread]
List: postgis-users
Subject: [postgis-users] ST_AsRaster and gdal_translate
From: Guillaume Drolet DRF <DROGU2 () intranet ! mrn ! gouv>
Date: 2012-07-26 15:11:56
Message-ID: 50115E3C.2070101 () intranet ! mrn ! gouv
[Download RAW message or body]
[Attachment #2 (multipart/alternative)]
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>
[Attachment #5 (text/html)]
<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=ISO-8859-1">
</head>
<body bgcolor="#FFFFFF" text="#000000">
<pre>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:
><i> If he does that he will not get the right value from his polygon layer...
</i>><i>
</i>><i> The right way to create an exportable raster from a polygon table it to:
</i>><i>
</i>><i> 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): \
</i>><i> </i>><i> CREATE TABLE rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m \
AS </i>><i> SELECT ST_AsRaster(hyd.geom, 1000.0, 1000.0, 0.0, 0.0, '32BUI'::text, \
nameOfTheColumnContainingTheValueToBurn ) AS rast </i>><i> FROM \
hydrology.bdtq_20k_hydro_so_gfsm_aea_subset AS hyd; </i>><i>
</i>><i> 2) Union the table.
</i>><i>
</i>><i> If your raster is very small you can do:
</i>><i>
</i>><i> CREATE TABLE rast2export AS
</i>><i> SELECT ST_Union(rast) rast
</i>><i> FROM rasters.bdtq_20k_hydro_so_gfsm_aea_subset
</i>><i>
</i>><i> If your raster is big, ST_Union is too slow. You have to use the method I \
describe in my blog </i>><i>
</i>><i> <a href="http://geospatialelucubrations.blogspot.ca/2012/07/a-slow-yet-100 \
0x-faster-alternative-to.html">http://geospatialelucubrations.blogspot.ca/2012/07/a-slow-yet-1000x-faster-alternative-to.html</a>
</i>><i>
</i>><i> 3) Then you can GDAL_translate using mode=1
</i>><i>
</i>><i> Sorry it is so complicated...
</i>><i>
</i>><i> 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. </i>><i>
</i>><i> Pierre
</i>><i>
</i>><i>
</i>><i>
</i>>><i> -----Original Message-----
</i>>><i> From: <a \
href="http://postgis.refractions.net/mailman/listinfo/postgis-users">postgis-users-bounces \
at postgis.refractions.net</a> [<a class="moz-txt-link-freetext" \
href="mailto:postgis-users">mailto:postgis-users</a>- </i>>><i> <a \
href="http://postgis.refractions.net/mailman/listinfo/postgis-users">bounces at \
postgis.refractions.net</a>] On Behalf Of Bborie Park </i>>><i> Sent: \
Wednesday, July 25, 2012 4:23 PM </i>>><i> To: <a \
href="http://postgis.refractions.net/mailman/listinfo/postgis-users">postgis-users at \
postgis.refractions.net</a> </i>>><i> Subject: Re: [postgis-users] ST_AsRaster \
and gdal_translate </i>>><i>
</i>>><i> For exporting, instead of using the rasters in
</i>>><i> "bdtq_20k_hydro_so_gfsm_aea_subset_1000m", you may want to try
</i>>><i> something
</i>>><i> like...
</i>>><i>
</i>>><i> CREATE TABLE foo AS
</i>>><i> SELECT 1 AS rid, ST_AsRaster((
</i>>><i> SELECT
</i>>><i> ST_Collect(geom)
</i>>><i> FROM hydrology.bdtq_20k_hydro_so_gfsm_aea_subset
</i>>><i> ), 1000.0, 1000.0
</i>>><i> ) AS rast
</i>>><i>
</i>>><i> The above creates one raster using the collection of geometries, \
which </i>>><i> can then be exported with gdal_translate.
</i>>><i>
</i>>><i> -bborie
</i>>><i>
</i>>><i> On 07/25/2012 12:57 PM, Guillaume Drolet DRF wrote:
</i>>>><i> Hi,
</i>>>><i>
</i>>>><i> I have been struggling with this since yesterday and now I think \
it is </i>>>><i> time for me to request some help from the community :
</i>>>><i>
</i>>>><i> I am trying to convert a geometry table (polygons) to a raster \
that will </i>>>><i> have 1000 m-pixels, and then export that raster to a \
TIF file using </i>>>><i> gdal_translate.
</i>>>><i> I've tried loads of different queries but nothing seems to \
work... Here </i>>>><i> are examples of some the queries I tried:
</i>>>><i>
</i>>>><i> CREATE TABLE rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m AS
</i>>>><i> SELECT ST_AsRaster(hyd.geom, 1000.0, 1000.0) AS rast
</i>>>><i> FROM hydrology.bdtq_20k_hydro_so_gfsm_aea_subset AS hyd;
</i>>>><i>
</i>>>><i> ALTER TABLE rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m ADD
</i>>><i> COLUMN
</i>>>><i> rid serial PRIMARY KEY;
</i>>>><i>
</i>>>><i> CREATE INDEX
</i>>>><i> ON rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m
</i>>>><i> USING gist(st_convexhull(rast));
</i>>>><i>
</i>>>><i> SELECT AddRasterConstraints('rasters'::name,
</i>>>><i> 'bdtq_20k_hydro_so_gfsm_aea_subset_1000m'::name,
</i>>>><i> 'rast'::name);
</i>>>><i>
</i>>>><i> Then at the command line:
</i>>>><i>
</i>>>>><i> gdal_translate "PG:host=localhost port=5432 dbname=testspdb
</i>>>><i> user=postgres password=*** schema=rasters
</i>>>><i> table=bdtq_20k_hydro_so_gfsm_aea_subset_1000m mode=2"
</i>>>><i> E:\temp\hydro1000m.tif
</i>>>><i>
</i>>>><i> I get this error:
</i>>>><i>
</i>>>><i> ERROR 1: Error, the table
</i>>>><i> rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m contains tiles \
with </i>>>><i> different size, and irregular blocking is not supported yet
</i>>>><i> GDALOpen failed - 1
</i>>>><i> Error, the table rasters.bdtq_20k_hydro_so_gfsm_aea_subset_1000m
</i>>>><i> contains tiles
</i>>>><i> with different size, and irregular blocking is not supported yet
</i>>>><i>
</i>>>><i> My problem seems to be related to the irregular blocking of my \
raster. I </i>>>><i> though I could create a raster with regular block and \
pixel size but I </i>>>><i> didn't figure how to do this yet. I tried to \
achieve this by creating an </i>>>><i> empty raster:
</i>>>><i>
</i>>>><i> CREATE TABLE rasters.test_rasterize (
</i>>>><i> rid serial PRIMARY KEY,
</i>>>><i> rast raster
</i>>>><i> );
</i>>>><i>
</i>>>><i> INSERT INTO rasters.test_rasterize VALUES (
</i>>>><i> 1,
</i>>>><i> ST_MakeEmptyRaster(45, 62, 1828552.00, 1422569.00, 1000.0, \
1000.0, </i>>>><i> 0.0, 0.0, 3175)
</i>>>><i> );
</i>>>><i> -- 45 and 62 are the approximate numbers of 1000 m pixels in x \
and y </i>>>><i> that fit inside the extent of my geometry layer
</i>>>><i> UPDATE rasters.test_rasterize
</i>>>><i> SET rast = ST_AsRaster(geom, 1000.0, 1000.0)
</i>>>><i> FROM hydrology.bdtq_20k_hydro_so_gfsm_aea_subset;
</i>>>><i>
</i>>>><i>
</i>>>><i> When I call gdal_translate with rasters.test_rasterize, I get \
this: </i>>>><i>
</i>>>><i> Input file size is 2, 4
</i>>>><i> 0...10...20...30...40...50...60...70...80...90...100 - done.
</i>>>><i>
</i>>>><i> I'm getting desperate and will greatly appreciate any hints!
</i>>>><i>
</i>>>><i> If that can help, I'm using PostGIS 2.0.0 r9605 and GDAL 1.9.0
</i>>>><i>
</i>>>><i> Many thanks,
</i>>>><i>
</i>>>><i> Guillaume
</i>>>><i>
</i>>>><i> _______________________________________________
</i>>>><i> postgis-users mailing list
</i>>>><i> <a \
href="http://postgis.refractions.net/mailman/listinfo/postgis-users">postgis-users at \
postgis.refractions.net</a> </i>>>><i> <a \
href="http://postgis.refractions.net/mailman/listinfo/postgis-users">http://postgis.refractions.net/mailman/listinfo/postgis-users</a>
</i>>>><i>
</i>>><i>
</i>>><i> --
</i>>><i> Bborie Park
</i>>><i> Programmer
</i>>><i> Center for Vectorborne Diseases
</i>>><i> UC Davis
</i>>><i> 530-752-8380
</i>>><i> <a href="http://postgis.refractions.net/mailman/listinfo/postgis-users">bkpark \
at ucdavis.edu</a> </i>>><i> _______________________________________________
</i>>><i> postgis-users mailing list
</i>>><i> <a href="http://postgis.refractions.net/mailman/listinfo/postgis-users">postgis-users \
at postgis.refractions.net</a> </i>>><i> <a \
href="http://postgis.refractions.net/mailman/listinfo/postgis-users">http://postgis.refractions.net/mailman/listinfo/postgis-users</a>
</i>><i> _______________________________________________
</i>><i> postgis-users mailing list
</i>><i> <a href="http://postgis.refractions.net/mailman/listinfo/postgis-users">postgis-users \
at postgis.refractions.net</a> </i>><i> <a \
href="http://postgis.refractions.net/mailman/listinfo/postgis-users">http://postgis.refractions.net/mailman/listinfo/postgis-users</a>
</i>><i>
</i>
--
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
<a href="http://postgis.refractions.net/mailman/listinfo/postgis-users">bkpark at \
ucdavis.edu</a></pre> </body>
</html>
_______________________________________________
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic