[postgis-users] ST_Buffer + grid problem

Pierre Racine Pierre.Racine at sbf.ulaval.ca
Thu Mar 22 09:09:14 PDT 2012


Hum... That works here. Could you do:

SELECT PostGIS_Full_Version();

> -----Original Message-----
> From: manohar.kaul at gmail.com [mailto:manohar.kaul at gmail.com] On Behalf
> Of Ed Linde
> Sent: Thursday, March 22, 2012 11:45 AM
> To: Pierre Racine
> Cc: PostGIS Users Discussion; Bborie Park (bkpark at ucdavis.edu)
> Subject: Re: [postgis-users] ST_Buffer + grid problem
> 
> 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
> 			>       >       >
> 			>       >       >
> 			>       >
> 			>       >
> 			>       >
> 			>
> 			>
> 			>
> 
> 
> 
> 
> 




More information about the postgis-users mailing list