[postgis-users] ST_Buffer + grid problem
Pierre Racine
Pierre.Racine at sbf.ulaval.ca
Thu Mar 22 09:59:02 PDT 2012
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
> > > > >
> > > > >
> > > >
> > > >
> > > >
> > >
> > >
> > >
> >
> >
> >
> >
> >
>
>
>
More information about the postgis-users
mailing list