[postgis-users] ST_Buffer + grid problem

Ed Linde edolinde at gmail.com
Thu Mar 22 10:19:35 PDT 2012


I can try that ... but I am away from the office now. So will have to give
it a shot tomorrow morning! :)
Will let you know if that fixes it. But should I download the source and
build GDAL 1.9?

On Thu, Mar 22, 2012 at 5:59 PM, Pierre Racine
<Pierre.Racine at sbf.ulaval.ca>wrote:

> What if you upgrade GDAL to 1.9?
>
> > -----Original Message-----
> > From: manohar.kaul at gmail.com [mailto:manohar.kaul at gmail.com] On Behalf
> > Of Ed Linde
> > Sent: Thursday, March 22, 2012 12:11 PM
> > To: Pierre Racine
> > Cc: PostGIS Users Discussion; Bborie Park (bkpark at ucdavis.edu)
> > Subject: Re: [postgis-users] ST_Buffer + grid problem
> >
> > "POSTGIS="2.0.0alpha7SVN" GEOS="3.3.0-CAPI-1.7.0" PROJ="Rel. 4.7.1, 23
> > September 2009" GDAL="GDAL 1.7.3, released 2010/11/10" LIBXML="2.7.8"
> > USE_STATS"
> > I used this because I remember I wanted GEOS/GDAL support and while
> installing
> > it there was some easy option of getting raster support. I think someone
> > suggested using Postgis 2.0!
> >
> >
> > On Thu, Mar 22, 2012 at 5:09 PM, Pierre Racine <
> Pierre.Racine at sbf.ulaval.ca>
> > wrote:
> >
> >
> >       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
> >       >                       >       >       >
> >       >                       >       >       >
> >       >                       >       >
> >       >                       >       >
> >       >                       >       >
> >       >                       >
> >       >                       >
> >       >                       >
> >       >
> >       >
> >       >
> >       >
> >       >
> >
> >
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20120322/8ae20559/attachment.html>


More information about the postgis-users mailing list