I can try that ... but I am away from the office now. So will have to give it a shot tomorrow morning! :)<div>Will let you know if that fixes it. But should I download the source and build GDAL 1.9?<br><br><div class="gmail_quote">
On Thu, Mar 22, 2012 at 5:59 PM, Pierre Racine <span dir="ltr"><<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
What if you upgrade GDAL to 1.9?<br>
<div class="im HOEnZb"><br>
> -----Original Message-----<br>
> From: <a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a> [mailto:<a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a>] On Behalf<br>
> Of Ed Linde<br>
</div><div class="HOEnZb"><div class="h5">> Sent: Thursday, March 22, 2012 12:11 PM<br>
> To: Pierre Racine<br>
> Cc: PostGIS Users Discussion; Bborie Park (<a href="mailto:bkpark@ucdavis.edu">bkpark@ucdavis.edu</a>)<br>
> Subject: Re: [postgis-users] ST_Buffer + grid problem<br>
><br>
> "POSTGIS="2.0.0alpha7SVN" GEOS="3.3.0-CAPI-1.7.0" PROJ="Rel. 4.7.1, 23<br>
> September 2009" GDAL="GDAL 1.7.3, released 2010/11/10" LIBXML="2.7.8"<br>
> USE_STATS"<br>
> I used this because I remember I wanted GEOS/GDAL support and while installing<br>
> it there was some easy option of getting raster support. I think someone<br>
> suggested using Postgis 2.0!<br>
><br>
><br>
> On Thu, Mar 22, 2012 at 5:09 PM, Pierre Racine <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>><br>
> wrote:<br>
><br>
><br>
>       Hum... That works here. Could you do:<br>
><br>
>       SELECT PostGIS_Full_Version();<br>
><br>
><br>
>       > -----Original Message-----<br>
>       > From: <a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a> [mailto:<a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a>]<br>
> On Behalf<br>
>       > Of Ed Linde<br>
><br>
>       > Sent: Thursday, March 22, 2012 11:45 AM<br>
>       > To: Pierre Racine<br>
><br>
>       > Cc: PostGIS Users Discussion; Bborie Park (<a href="mailto:bkpark@ucdavis.edu">bkpark@ucdavis.edu</a>)<br>
>       > Subject: Re: [postgis-users] ST_Buffer + grid problem<br>
>       ><br>
>       > Hi Pierre,<br>
>       > Here goes the error reproduced (hope it helps): If this cannot be<br>
> resolved is<br>
>       > there another work around? I am a bit under the pump here to get this<br>
> grid done<br>
>       > :(<br>
>       > So the SRID = 4326 and this is the polygon/bounding box for entire<br>
> road network<br>
>       > of OSM.<br>
>       ><br>
>       > select ST_AsRaster(<br>
> ST_GeomFromText('POLYGON((8.07734039737749<br>
>       > 54.4984986588244,8.07734039737749<br>
> 57.7505109647578,15.1919565742587<br>
>       > 57.7505109647578,15.1919565742587<br>
> 54.4984986588244,8.07734039737749<br>
>       > 54.4984986588244))'), 0.000036, 0.000036);<br>
>       ><br>
>       > ERROR:  rt_raster_gdal_rasterize: Unable to add band to GDALDataset<br>
>       ><br>
>       > ********** Error **********<br>
>       ><br>
>       > ERROR: rt_raster_gdal_rasterize: Unable to add band to GDALDataset<br>
>       > SQL state: XX000<br>
>       ><br>
>       ><br>
>       ><br>
>       ><br>
>       > On Thu, Mar 22, 2012 at 4:00 PM, Ed Linde <<a href="mailto:edolinde@gmail.com">edolinde@gmail.com</a>><br>
> wrote:<br>
>       ><br>
>       ><br>
>       >       I tried a few roads at random and it doesn't reproduce that error,<br>
>       > though I am a bit suspicious about the actual<br>
>       >       query. This part in particular<br>
>       ><br>
>       >       (SELECT ST_AsRaster(ST_Extent(way_geom)::geometry, 0.000036,<br>
>       >       > 0.000036) rast<br>
>       >       >                          FROM buffers<br>
>       >       >                         ) foo1<br>
>       ><br>
>       >       Are we saying that we want to compute the extent of each and<br>
> every<br>
>       > geometry in my table? Or should we not<br>
>       >       be getting the maximum extents for ALL the geometries?<br>
>       >       So when I say something like "FROM buffers where osm_id in (...<br>
> list of<br>
>       > 10 IDs)... " this query runs forever!<br>
>       >       Eventhough checking them one by one returns with results<br>
> immediately.<br>
>       >       So I am not entirely sure how to do this or how I can reproduce this<br>
> bug<br>
>       > to help you guys.<br>
>       >       Any suggestions?<br>
>       ><br>
>       ><br>
>       >       On Thu, Mar 22, 2012 at 3:23 PM, Ed Linde <<a href="mailto:edolinde@gmail.com">edolinde@gmail.com</a>><br>
>       > wrote:<br>
>       ><br>
>       ><br>
>       >               Yeah only problem being that I don't know which geometry is<br>
>       > the one that fails... its a table with 300K geometries.<br>
>       >               Have never written any PgSQL, so this I guess would require a<br>
>       > loop to check which OSM ID its failing on. Unless<br>
>       >               you got a good idea to pin-point this error?<br>
>       ><br>
>       ><br>
>       >               On Thu, Mar 22, 2012 at 3:16 PM, Pierre Racine<br>
>       > <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>> wrote:<br>
>       ><br>
>       ><br>
>       >                       What I mean is to replace:<br>
>       ><br>
>       ><br>
>       >                       SELECT ST_AsRaster(ST_Extent(way_geom)::geometry,<br>
>       > 0.000036, 0.000036) rast FROM buffers  where osm_id = 94695311<br>
>       ><br>
>       ><br>
>       >                       by something like:<br>
>       ><br>
>       >                       SELECT<br>
>       > ST_AsRaster(ST_GeomFromText('POLYGON((11.1539754732244<br>
>       > 55.5772503545236,11.1540114732244<br>
> 55.5772503545236,11.1540114732244<br>
>       > 55.5772863545236,11.1539754732244<br>
> 55.5772863545236,11.1539754732244<br>
>       > 55.5772503545236))'), 0.000036, 0.000036)<br>
>       ><br>
>       >                       with one of the geometry reproducing the bug so we<br>
>       > don't need your dataset to test and reproduce.<br>
>       ><br>
>       ><br>
>       >                       Pierre<br>
>       ><br>
>       >                       > -----Original Message-----<br>
>       >                       > From: <a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a><br>
>       > [mailto:<a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a>] On Behalf<br>
>       >                       > Of Ed Linde<br>
>       ><br>
>       >                       > Sent: Thursday, March 22, 2012 9:59 AM<br>
>       >                       > To: Pierre Racine<br>
>       ><br>
>       >                       > Cc: PostGIS Users Discussion; Bborie Park<br>
>       > (<a href="mailto:bkpark@ucdavis.edu">bkpark@ucdavis.edu</a>)<br>
>       >                       > Subject: Re: [postgis-users] ST_Buffer + grid problem<br>
>       >                       ><br>
>       >                       > sure. Did you mean ST_GeomfromText or show it as<br>
>       > text?<br>
>       >                       ><br>
>       >                       > SELECT ST_GeomfromText ((gvxy).geom), ((gvxy).x -<br>
>       > 1) * rwidth + (gvxy).y gridid<br>
>       >                       > FROM (SELECT ST_PixelAsPolygons(rast) gvxy,<br>
>       > ST_Width(rast) rwidth<br>
>       >                       >             FROM (SELECT<br>
>       > ST_AsRaster(ST_Extent(way_geom)::geometry, 0.000036,<br>
>       >                       > 0.000036) rast<br>
>       >                       >                          FROM buffers<br>
>       >                       >                          where osm_id = 94695311<br>
>       >                       >                         ) foo1<br>
>       >                       >            ) foo2;<br>
>       >                       ><br>
>       >                       > ERROR:  parse error - invalid geometry<br>
>       >                       > HINT:  "010300000001000000050000005" <-- parse<br>
>       > error at position 27 within<br>
>       >                       > geometry<br>
>       >                       ><br>
>       >                       > ********** Error **********<br>
>       >                       ><br>
>       >                       > ERROR: parse error - invalid geometry<br>
>       >                       > SQL state: XX000<br>
>       >                       > Hint: "010300000001000000050000005" <-- parse<br>
>       > error at position 27 within<br>
>       >                       > geometry<br>
>       >                       ><br>
>       >                       ><br>
>       >                       ><br>
>       >                       > On Thu, Mar 22, 2012 at 2:53 PM, Pierre Racine<br>
>       > <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>><br>
>       >                       > wrote:<br>
>       >                       ><br>
>       >                       ><br>
>       >                       >       I guess this is related to a GDAL problem where<br>
>       > you cannot create a<br>
>       >                       > raster using ST_AsRaster with pixel size smaller than<br>
>       > 1. Was this fixed Bborie?<br>
>       >                       ><br>
>       >                       >       Ed could you reduce the list of geometries passed<br>
>       > to ST_Extent to one<br>
>       >                       > and write it as ST_GeomfromText so we can debug<br>
>       > this more easily?<br>
>       >                       ><br>
>       >                       ><br>
>       >                       >       Pierre<br>
>       >                       ><br>
>       >                       >       > -----Original Message-----<br>
>       >                       >       > From: <a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a><br>
>       > [mailto:<a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a>]<br>
>       >                       > On Behalf<br>
>       >                       >       > Of Ed Linde<br>
>       >                       ><br>
>       >                       >       > Sent: Thursday, March 22, 2012 9:45 AM<br>
>       >                       >       > To: Pierre Racine<br>
>       >                       >       > Cc: PostGIS Users Discussion<br>
>       >                       >       > Subject: Re: [postgis-users] ST_Buffer + grid<br>
>       > problem<br>
>       >                       >       ><br>
>       >                       >       > Hi Pierre,<br>
>       >                       ><br>
>       >                       >       > So I am trying to have a 4meter by 4meter cell<br>
>       > and I converted<br>
>       >                       > 4meters to<br>
>       >                       >       > degrees, because I think ST_Extent will compute<br>
>       > it in degrees as my<br>
>       >                       > geometry is<br>
>       >                       >       > in srid 4326. But I get this GDAL error which I<br>
>       > have no clue what it<br>
>       >                       > means. Any<br>
>       >                       >       > ideas how I can fix it? Also I noticed that if I set it<br>
>       > to 1.0, 1.0 I get<br>
>       >                       >       > 21 rows but the gridid isn't continuous either<br>
>       > (which is weird).<br>
>       >                       >       ><br>
>       >                       >       > CREATE TABLE vectorgrid AS<br>
>       >                       >       > SELECT (gvxy).geom, ((gvxy).x - 1) * rwidth +<br>
>       > (gvxy).y gridid FROM<br>
>       >                       > (SELECT<br>
>       >                       >       > ST_PixelAsPolygons(rast) gvxy, ST_Width(rast)<br>
>       > rwidth<br>
>       >                       >       >             FROM (SELECT<br>
>       > ST_AsRaster(ST_Extent(way_geom)::geometry,<br>
>       >                       > 0.000036,<br>
>       >                       >       > 0.000036) rast<br>
>       >                       >       >                          FROM buffers<br>
>       >                       >       >                         ) foo1<br>
>       >                       >       >            ) foo2;<br>
>       >                       >       ><br>
>       >                       >       > ERROR:  rt_raster_gdal_rasterize: Unable to add<br>
>       > band to GDALDataset<br>
>       >                       >       ><br>
>       >                       >       > ********** Error **********<br>
>       >                       >       ><br>
>       >                       >       > ERROR: rt_raster_gdal_rasterize: Unable to add<br>
>       > band to GDALDataset<br>
>       >                       > SQL state:<br>
>       >                       >       > XX000<br>
>       >                       >       ><br>
>       >                       >       ><br>
>       >                       >       ><br>
>       >                       >       ><br>
>       >                       >       > On Thu, Mar 22, 2012 at 2:13 PM, Pierre Racine<br>
>       >                       > <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>><br>
>       >                       >       > wrote:<br>
>       >                       >       ><br>
>       >                       >       ><br>
>       >                       >       >       You can control how the grid is aligned by<br>
>       > using more ST_AsRaster<br>
>       >                       >       > parameters. See:<br>
>       >                       >       ><br>
>       >                       >       ><br>
>       > <a href="http://postgis.refractions.net/documentation/manual-" target="_blank">http://postgis.refractions.net/documentation/manual-</a><br>
>       >                       >       > svn/RT_ST_AsRaster.html<br>
>       >                       >       ><br>
>       >                       >       >       If you want it to align on your point, just align<br>
>       > it on your points...<br>
>       >                       >       ><br>
>       >                       >       >       Pierre<br>
>       >                       >       ><br>
>       >                       >       ><br>
>       >                       >       >       > -----Original Message-----<br>
>       >                       >       >       > From: <a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a><br>
>       >                       > [mailto:<a href="mailto:manohar.kaul@gmail.com">manohar.kaul@gmail.com</a>]<br>
>       >                       >       > On Behalf<br>
>       >                       >       >       > Of Ed Linde<br>
>       >                       >       >       > Sent: Thursday, March 22, 2012 9:05 AM<br>
>       >                       >       >       > To: PostGIS Users Discussion<br>
>       >                       >       >       > Cc: Pierre Racine<br>
>       >                       >       >       > Subject: Re: [postgis-users] ST_Buffer + grid<br>
>       > problem<br>
>       >                       >       >       ><br>
>       >                       >       >       > Hi Pierre,<br>
>       >                       >       >       > Thanks for the grid idea last night. I have a<br>
>       > problem though, I<br>
>       >                       > want to<br>
>       >                       >       > have a 4m<br>
>       >                       >       >       > by 4m cell sized grid over my buffer<br>
>       > geometries which are in SRID<br>
>       >                       > =<br>
>       >                       >       > 4326. And<br>
>       >                       >       >       > then I have another data set of 2D points<br>
>       > which are also in SRID<br>
>       >                       >       > 4326...<br>
>       >                       >       >       > wondering how I can figure out from just a<br>
>       > given (lat,long) which<br>
>       >                       > cell<br>
>       >                       >       > ID in the<br>
>       >                       >       >       > buffer grid it would belong to? I am thinking<br>
>       > that I might run into<br>
>       >                       >       > "alignment<br>
>       >                       >       >       > issues" because shouldn't the extent of the<br>
>       > grid be exactly the<br>
>       >                       > same<br>
>       >                       >       > on the 2D<br>
>       >                       >       >       > point cloud as well, so that the grid cell IDs<br>
>       > will match? Unless I<br>
>       >                       > am<br>
>       >                       >       > missing<br>
>       >                       >       >       > something here?<br>
>       >                       >       >       ><br>
>       >                       >       >       > Cheers,<br>
>       >                       >       >       > Ed<br>
>       >                       >       >       ><br>
>       >                       >       >       ><br>
>       >                       >       >       ><br>
>       >                       >       >       > On Wed, Mar 21, 2012 at 8:46 PM, Pierre<br>
>       > Racine<br>
>       >                       >       > <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>><br>
>       >                       >       >       > wrote:<br>
>       >                       >       >       ><br>
>       >                       >       >       ><br>
>       >                       >       >       >       > I am not sure if this is possible, but I<br>
>       > have computed (using<br>
>       >                       >       > ST_Buffer)<br>
>       >                       >       >       > a sort of<br>
>       >                       >       >       >       > buffer around several LINESTRINGs.<br>
>       > Now I would like to lay<br>
>       >                       > some<br>
>       >                       >       > sort<br>
>       >                       >       >       > of grid of<br>
>       >                       >       >       >       > say 1m^2 cell size on top of this<br>
>       > collection of geometries<br>
>       >                       > and<br>
>       >                       >       > then<br>
>       >                       >       >       > compute the<br>
>       >                       >       >       >       > intersection... so in the end I would<br>
>       > like to for example,<br>
>       >                       > know<br>
>       >                       >       > that grid<br>
>       >                       >       >       > cell ID = 1<br>
>       >                       >       >       >       > intersects with buffers 1 and 10, grid<br>
>       > cell 2 intersects with<br>
>       >                       > buffers<br>
>       >                       >       > 7, 10<br>
>       >                       >       >       > etc..<br>
>       >                       >       >       >       ><br>
>       >                       >       >       >       > Is there a simple way of doing this in<br>
>       > postgis? Maybe<br>
>       >                       > someone<br>
>       >                       >       > could<br>
>       >                       >       >       > point me to<br>
>       >                       >       >       >       > some documentation of how I can<br>
>       > generate such a grid in<br>
>       >                       > postgis<br>
>       >                       >       > and<br>
>       >                       >       >       > maybe<br>
>       >                       >       >       >       > then I can use just ST_Intersect once I<br>
>       > have these two<br>
>       >                       >       > geometries?<br>
>       >                       >       >       ><br>
>       >                       >       >       ><br>
>       >                       >       >       >       With the raster type you can now easily<br>
>       > create a vector grid<br>
>       >                       > like<br>
>       >                       >       > this:<br>
>       >                       >       >       ><br>
>       >                       >       >       >       CREATE TABLE vectorgrid AS<br>
>       >                       >       >       >       SELECT (gvxy).geom, ((gvxy).x - 1) *<br>
>       > rwidth + (gvxy).y gridid<br>
>       >                       >       >       >       FROM (SELECT ST_PixelAsPolygons(rast)<br>
>       > gvxy, ST_Width(rast)<br>
>       >                       >       > rwidth<br>
>       >                       >       >       >                   FROM (SELECT<br>
>       >                       > ST_AsRaster(ST_Extent(geom)::geometry,<br>
>       >                       >       > 1.0,<br>
>       >                       >       >       > 1.0) rast<br>
>       >                       >       >       >                                FROM yourbuffertable<br>
>       >                       >       >       >                               ) foo1<br>
>       >                       >       >       >                  ) foo2;<br>
>       >                       >       >       ><br>
>       >                       >       >       >       Make sure a spatial index exist on both<br>
>       > tables:<br>
>       >                       >       >       ><br>
>       >                       >       >       >       CREATE INDEX<br>
>       > yourbuffertable_geom_idx  ON<br>
>       >                       > yourbuffertable<br>
>       >                       >       > USING<br>
>       >                       >       >       > gist  (geom);<br>
>       >                       >       >       >       CREATE INDEX vectorgrid _geom_idx<br>
>       > ON vectorgrid USING<br>
>       >                       > gist<br>
>       >                       >       >       > (geom);<br>
>       >                       >       >       ><br>
>       >                       >       >       >       You can then perform a normal intersect<br>
>       > query:<br>
>       >                       >       >       ><br>
>       >                       >       >       >       CREATE TABLE interresult AS<br>
>       >                       >       >       >       SELECT b.bufferid, g.gridid,<br>
>       > ST_Intersection(g.geom, b.geom)<br>
>       >                       >       > geom<br>
>       >                       >       >       >       FROM vectorgrid g, yourbuffertable b<br>
>       >                       >       >       >       WHERE ST_Intersects(g.geom, b.geom);<br>
>       >                       >       >       ><br>
>       >                       >       >       >       Pierre<br>
>       >                       >       >       ><br>
>       > _______________________________________________<br>
>       >                       >       >       >       postgis-users mailing list<br>
>       >                       >       >       >       <a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br>
>       >                       >       >       ><br>
>       > <a href="http://postgis.refractions.net/mailman/listinfo/postgis-users" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a><br>
>       >                       >       >       ><br>
>       >                       >       >       ><br>
>       >                       >       ><br>
>       >                       >       ><br>
>       >                       >       ><br>
>       >                       ><br>
>       >                       ><br>
>       >                       ><br>
>       ><br>
>       ><br>
>       ><br>
>       ><br>
>       ><br>
><br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div>