Rasterize a vector

[Paolo Cavallini]

Hi all.
I'm testing the ST_AsRaster command to convert a vector of adjacent polygons into a
raster with different cell values according to one of the variables. I want the
result to align to an existing raster. Sounds an easy task, but it is proving elusive.

create table rasterizzo6
 as select id as rid,
 St_asRaster(geom, 100, 100, '4BUI')
 as rast from province;

creates a raster with irregular tiles, so it cannot be read by QGIS, and it is not
aligned to an existing rasters.

Are there more detailed instructions available?
Thanks.
--
Paolo Cavallini - Faunalia
www.faunalia.eu
Full contact details at www.faunalia.eu/pc
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Pierre Racine-2]

There are more parameters to ST_AsRaster to align the raster you produce. You can align on an existing raster or provide all the parameters (scalex, scaley, gridx, gridy) to align it.

See http://postgis.refractions.net/documentation/manual-svn/RT_ST_AsRaster.html

Pierre

> -----Original Message-----
> From: [hidden email] [mailto:postgis-users-
> [hidden email]] On Behalf Of Paolo Cavallini
> Sent: Tuesday, May 29, 2012 9:42 AM
> To: PostGIS Users Discussion
> Subject: [postgis-users] Rasterize a vector
>
> Hi all.
> I'm testing the ST_AsRaster command to convert a vector of adjacent polygons
> into a
> raster with different cell values according to one of the variables. I want the
> result to align to an existing raster. Sounds an easy task, but it is proving elusive.
>
> create table rasterizzo6
>  as select id as rid,
>  St_asRaster(geom, 100, 100, '4BUI')
>  as rast from province;
>
> creates a raster with irregular tiles, so it cannot be read by QGIS, and it is not
> aligned to an existing rasters.
>
> Are there more detailed instructions available?
> Thanks.
> --
> Paolo Cavallini - Faunalia
> www.faunalia.eu
> Full contact details at www.faunalia.eu/pc
> _______________________________________________
> postgis-users mailing list
> [hidden email]
> http://postgis.refractions.net/mailman/listinfo/postgis-users

_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Paolo Cavallini]

Il 29/05/2012 16:12, Pierre Racine ha scritto:
> There are more parameters to ST_AsRaster to align the raster you produce. You can align on an existing raster or provide all the parameters (scalex, scaley, gridx, gridy) to align it.
>
> See http://postgis.refractions.net/documentation/manual-svn/RT_ST_AsRaster.html

yes, I have read it, but I think it is (should be) much easier to use an existing raster:
ST_AsRaster(geometry geom, raster ref...
thanks.
--
Paolo Cavallini - Faunalia
www.faunalia.eu
Full contact details at www.faunalia.eu/pc
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Bborie Park]

On 05/29/2012 07:20 AM, Paolo Cavallini wrote:
> Il 29/05/2012 16:12, Pierre Racine ha scritto:
>> There are more parameters to ST_AsRaster to align the raster you produce. You can align on an existing raster or provide all the parameters (scalex, scaley, gridx, gridy) to align it.
>>
>> See http://postgis.refractions.net/documentation/manual-svn/RT_ST_AsRaster.html
>
> yes, I have read it, but I think it is (should be) much easier to use an existing raster:
> ST_AsRaster(geometry geom, raster ref...
> thanks.

Read it again.  The very first variant listed on that page has the
ST_AsRaster(geom geometry, rast raster, ...)

-bborie

--
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
[hidden email]
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Paolo Cavallini]

Il 29/05/2012 19:03, Bborie Park ha scritto:
> Read it again.  The very first variant listed on that page has the
> ST_AsRaster(geom geometry, rast raster, ...)

Yes, that's what I was trying, without success. Let's say I have two adjacent
polygons, and I want a raster where the cells have the value of one of the fields of
the the polygons (say, the id). Moreover, I want the resulting raster to be aligned
to an existing raster (this is a very common task in GIS, easily done in GRASS, GDAL,
SAGA, etc.). I do not see an option in:

ST_AsRaster(geometry geom, raster ref, text pixeltype, double precision value=1,
double precision nodataval=0, boolean touched=false);

to write in each cell the value of the raster. Moreover, a simple

 CREATE TABLE rasterized as
   select b.rid, ST_AsRaster(a.geom, b.rast, '4BUI') as rast
   from province a, hdr b;

results in (N geometries * N tiles) records, whereas I needed (N tiles) records, and
for each tile a raster with a different value for each province.

All the best, and thanks.
--
Paolo Cavallini - Faunalia
www.faunalia.eu
Full contact details at www.faunalia.eu/pc
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Pierre Racine-2]

> Yes, that's what I was trying, without success. Let's say I have two adjacent
> polygons, and I want a raster where the cells have the value of one of the fields
> of
> the the polygons (say, the id). Moreover, I want the resulting raster to be aligned
> to an existing raster (this is a very common task in GIS, easily done in GRASS,
> GDAL,
> SAGA, etc.). I do not see an option in:
>
> ST_AsRaster(geometry geom, raster ref, text pixeltype, double precision
> value=1,
> double precision nodataval=0, boolean touched=false);
>
> to write in each cell the value of the raster. Moreover, a simple
>
>  CREATE TABLE rasterized as
>    select b.rid, ST_AsRaster(a.geom, b.rast, '4BUI') as rast
>    from province a, hdr b;
>
> results in (N geometries * N tiles) records, whereas I needed (N tiles) records,
> and
> for each tile a raster with a different value for each province.

There should not be N tiles. Just N rasters...

The most important is: At the end do you want all geometries to be burned into the same raster or you want one raster per geometries?

Pierre

_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Bborie Park]

On 05/29/2012 10:56 AM, Paolo Cavallini wrote:

> Il 29/05/2012 19:03, Bborie Park ha scritto:
>> Read it again.  The very first variant listed on that page has the
>> ST_AsRaster(geom geometry, rast raster, ...)
>
> Yes, that's what I was trying, without success. Let's say I have two adjacent
> polygons, and I want a raster where the cells have the value of one of the fields of
> the the polygons (say, the id). Moreover, I want the resulting raster to be aligned
> to an existing raster (this is a very common task in GIS, easily done in GRASS, GDAL,
> SAGA, etc.). I do not see an option in:
>
> ST_AsRaster(geometry geom, raster ref, text pixeltype, double precision value=1,
> double precision nodataval=0, boolean touched=false);
>
> to write in each cell the value of the raster. Moreover, a simple
>
>  CREATE TABLE rasterized as
>    select b.rid, ST_AsRaster(a.geom, b.rast, '4BUI') as rast
>    from province a, hdr b;
>
> results in (N geometries * N tiles) records, whereas I needed (N tiles) records, and
> for each tile a raster with a different value for each province.
>
> All the best, and thanks.

WITH ref AS (
        SELECT ST_MakeEmptyRaster(1, 1, 0, 0, 1, -1, 0, 0, 0) AS rast
)
SELECT
        ST_AsRaster(g.geom, ref.rast, '32BUI', g.id, 0)
FROM mygeomtable g
CROSS JOIN ref

That would convert your geometries into rasters, all of which have the
same alignment with pixel values set to the particular geometry's id.
So, for 20 input geometries, there would be 20 output rasters.

-bborie

--
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
[hidden email]
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Paolo Cavallini]

Il 29/05/2012 20:05, Pierre Racine ha scritto:

> The most important is: At the end do you want all geometries to be burned into the same raster or you want one raster per geometries?

same raster - see gdal_rasterize for an example.
thanks a lot.

--
Paolo Cavallini - Faunalia
www.faunalia.eu
Full contact details at www.faunalia.eu/pc
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Bborie Park]

On 05/29/2012 11:07 AM, Paolo Cavallini wrote:
> Il 29/05/2012 20:05, Pierre Racine ha scritto:
>
>> The most important is: At the end do you want all geometries to be burned into the same raster or you want one raster per geometries?
>
> same raster - see gdal_rasterize for an example.
> thanks a lot.
>

Run an ST_Union upon the set of output rasters.

-bborie

--
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
[hidden email]
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Paolo Cavallini]

Il 29/05/2012 20:06, Bborie Park ha scritto:

> WITH ref AS (
> SELECT ST_MakeEmptyRaster(1, 1, 0, 0, 1, -1, 0, 0, 0) AS rast
> )
> SELECT
> ST_AsRaster(g.geom, ref.rast, '32BUI', g.id, 0)
> FROM mygeomtable g
> CROSS JOIN ref
>
> That would convert your geometries into rasters, all of which have the
> same alignment with pixel values set to the particular geometry's id.
> So, for 20 input geometries, there would be 20 output rasters.

but:

WITH ref AS (
SELECT ST_MakeEmptyRaster(1, 1, 0, 0, 1, -1, 0, 0, 3003) AS rast
)
SELECT
ST_AsRaster(g.geom, ref.rast, '32BUI', g.id, 0)
FROM province g
CROSS JOIN ref;
ERROR:  rt_raster_gdal_rasterize: Unable to add band to GDALDataset
CONTEXT: PL/pgSQL function "st_asraster" line 26 at RETURN

All the best.
--
Paolo Cavallini - Faunalia
www.faunalia.eu
Full contact details at www.faunalia.eu/pc
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Pierre Racine-2]

Basically (I use the old sub query style):

SELECT ST_Union(ST_AsRaster(g.geom, ref.rast, '32BUI', g.id, 0))) rast
FROM (SELECT ST_MakeEmptyRaster(1, 1, 0, 0, 1, -1, 0, 0, 0) AS rast) ref, mygeomtable g

however if your raster is high resolution ST_Union might be impractically slow.

The other technique untill we get a faster ST_Union is to use ST_MapAlgebraFct and write your own function returning the geometry value from the intersection with the centroid of the pixel. See this example http://postgis.refractions.net/pipermail/postgis-users/2012-April/033875.html

I agree this is relatively cumbersome but it should be much faster than ST_Union.

Pierre

> -----Original Message-----
> From: [hidden email] [mailto:postgis-users-
> [hidden email]] On Behalf Of Bborie Park
> Sent: Tuesday, May 29, 2012 2:11 PM
> To: [hidden email]
> Subject: Re: [postgis-users] Rasterize a vector
>
> On 05/29/2012 11:07 AM, Paolo Cavallini wrote:
> > Il 29/05/2012 20:05, Pierre Racine ha scritto:
> >
> >> The most important is: At the end do you want all geometries to be burned
> into the same raster or you want one raster per geometries?
> >
> > same raster - see gdal_rasterize for an example.
> > thanks a lot.
> >
>
> Run an ST_Union upon the set of output rasters.
>
> -bborie
>
> --
> Bborie Park
> Programmer
> Center for Vectorborne Diseases
> UC Davis
> 530-752-8380
> [hidden email]
> _______________________________________________
> postgis-users mailing list
> [hidden email]
> http://postgis.refractions.net/mailman/listinfo/postgis-users

_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Bborie Park]

On 05/29/2012 11:20 AM, Paolo Cavallini wrote:

> Il 29/05/2012 20:06, Bborie Park ha scritto:
>
>> WITH ref AS (
>> SELECT ST_MakeEmptyRaster(1, 1, 0, 0, 1, -1, 0, 0, 0) AS rast
>> )
>> SELECT
>> ST_AsRaster(g.geom, ref.rast, '32BUI', g.id, 0)
>> FROM mygeomtable g
>> CROSS JOIN ref
>>
>> That would convert your geometries into rasters, all of which have the
>> same alignment with pixel values set to the particular geometry's id.
>> So, for 20 input geometries, there would be 20 output rasters.
>
> but:
>
> WITH ref AS (
> SELECT ST_MakeEmptyRaster(1, 1, 0, 0, 1, -1, 0, 0, 3003) AS rast
> )
> SELECT
> ST_AsRaster(g.geom, ref.rast, '32BUI', g.id, 0)
> FROM province g
> CROSS JOIN ref;
> ERROR:  rt_raster_gdal_rasterize: Unable to add band to GDALDataset
> CONTEXT: PL/pgSQL function "st_asraster" line 26 at RETURN
>

Strange.  That works for me on a variety of my geometry tables...

WITH ref AS (
        SELECT ST_MakeEmptyRaster(1, 1, 0, 0, 1, -1, 0, 0, 4326) AS rast
)
SELECT
        (ST_SummaryStats(ST_AsRaster(s.shape, ref.rast, '32BUI', s.id, -1))).*
FROM site s
CROSS JOIN ref
LIMIT 10

Can you provide the output of

SELECT postgis_full_version();

-bborie
--
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
[hidden email]
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Paolo Cavallini]

Il 29/05/2012 20:25, Pierre Racine ha scritto:
> Basically (I use the old sub query style):
>
> SELECT ST_Union(ST_AsRaster(g.geom, ref.rast, '32BUI', g.id, 0))) rast
> FROM (SELECT ST_MakeEmptyRaster(1, 1, 0, 0, 1, -1, 0, 0, 0) AS rast) ref, mygeomtable g

same error: ERROR:  rt_raster_gdal_rasterize: Unable to add band to GDALDataset
CONTEXT: PL/pgSQL function "st_asraster" line 26 at RETURN

It's becoming obvious that it is much more practical and faster to do the
rasterization outside the DB, at lest for now.
Thanks a lot for the insights.
All the best.
--
Paolo Cavallini - Faunalia
www.faunalia.eu
Full contact details at www.faunalia.eu/pc
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Paolo Cavallini]

Il 29/05/2012 20:33, Bborie Park ha scritto:

> Can you provide the output of
>
> SELECT postgis_full_version();
>
> -bborie

Sure:

 POSTGIS="2.0.0 r9605" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23 September 2009"
GDAL="GDAL 1.9.0, released 2011/12/29" LIBXML="2.7.8" RASTER

All the best, and many thanks.
--
Paolo Cavallini - Faunalia
www.faunalia.eu
Full contact details at www.faunalia.eu/pc
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Pierre Racine-2]

> same error: ERROR:  rt_raster_gdal_rasterize: Unable to add band to
> GDALDataset
> CONTEXT: PL/pgSQL function "st_asraster" line 26 at RETURN

Isn't that because he has to set GDAL_DATA or PROJSO?

Pierre
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Bborie Park]

On 05/29/2012 11:35 AM, Paolo Cavallini wrote:

> Il 29/05/2012 20:33, Bborie Park ha scritto:
>
>> Can you provide the output of
>>
>> SELECT postgis_full_version();
>>
>> -bborie
>
> Sure:
>
>  POSTGIS="2.0.0 r9605" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23 September 2009"
> GDAL="GDAL 1.9.0, released 2011/12/29" LIBXML="2.7.8" RASTER
>
> All the best, and many thanks.

Strange indeed.  The error itself arises from adding an empty band to a
GDAL MEM dataset, which could be caused by a memory issue (lack of
system memory) or invalid pixel type.  But  your pixel type is correct...

-bborie
--
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
[hidden email]
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Sandro Santilli]

On Tue, May 29, 2012 at 11:58:02AM -0700, Bborie Park wrote:

> On 05/29/2012 11:35 AM, Paolo Cavallini wrote:
> > Il 29/05/2012 20:33, Bborie Park ha scritto:
> >
> >> Can you provide the output of
> >>
> >> SELECT postgis_full_version();
> >>
> >> -bborie
> >
> > Sure:
> >
> >  POSTGIS="2.0.0 r9605" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23 September 2009"
> > GDAL="GDAL 1.9.0, released 2011/12/29" LIBXML="2.7.8" RASTER
> >
> > All the best, and many thanks.
>
> Strange indeed.  The error itself arises from adding an empty band to a
> GDAL MEM dataset, which could be caused by a memory issue (lack of
> system memory) or invalid pixel type.  But  your pixel type is correct...

Sign of a memory leak ?

--strk;

  ,------o-.
  |   __/  |    Delivering high quality PostGIS 2.0 !
  |  / 2.0 |    http://strk.keybit.net - http://vizzuality.com
  `-o------'

_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Bborie Park]

On 05/29/2012 12:20 PM, Sandro Santilli wrote:

> On Tue, May 29, 2012 at 11:58:02AM -0700, Bborie Park wrote:
>> On 05/29/2012 11:35 AM, Paolo Cavallini wrote:
>>> Il 29/05/2012 20:33, Bborie Park ha scritto:
>>>
>>>> Can you provide the output of
>>>>
>>>> SELECT postgis_full_version();
>>>>
>>>> -bborie
>>>
>>> Sure:
>>>
>>>  POSTGIS="2.0.0 r9605" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23 September 2009"
>>> GDAL="GDAL 1.9.0, released 2011/12/29" LIBXML="2.7.8" RASTER
>>>
>>> All the best, and many thanks.
>>
>> Strange indeed.  The error itself arises from adding an empty band to a
>> GDAL MEM dataset, which could be caused by a memory issue (lack of
>> system memory) or invalid pixel type.  But  your pixel type is correct...
>
> Sign of a memory leak ?
>

Possible.  I have to run valgrind on one of my monster geometry tables
and see.  I always run valgrind on the regression tests though before I
commit but something could have slipped through.

-bborie

--
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
[hidden email]
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Paolo Cavallini]

Il 29/05/2012 20:37, Pierre Racine ha scritto:

> Isn't that because he has to set GDAL_DATA or PROJSO?

You mean this?
http://trac.osgeo.org/gdal/wiki/FAQInstallationAndBuilding#HowtosetGDAL_DATAvariable
Gdal is installed from deb, so it's in the path.
All the best, and thanks.
--
Paolo Cavallini - Faunalia
www.faunalia.eu
Full contact details at www.faunalia.eu/pc
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Neil Best]

Paolo, were you able to get something working satisfactorily?  If I read the thread correctly you wanted a separate raster output for each geometry in the input.  I am more interested in the case of creating a raster coverage for a set of geometries, US counties for example, using an ID as the value.  It sounds like this distinction is not important, rather there was a problem with the construction of any output from these two cases that prevented it from being read by GDAL.  Do I have that right?  What is the state of the art on this front?  I am reluctant to try it if it is going to involve repeated impact of my forehead and the wall.  Thanks.

Neil