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