Thanks I will try and see how I go with all your suggestions tomorrow... hopefully you will be around to help. :) I have to play around <div>more with ST_MapalgebraExpr to understand how what you say works... like what is a "single band" and "double band" raster? </div>
<div>Does that just mean that in a single band the raster cells will have a 1D value as opposed to a 2D (x,y) in the double band?</div><div>Once the intersection between the road buffer geometry and this new raster is done, then the point cloud stuff I will have to write some code in C/C++ to handle that, so that shouldn't be a problem. One last thing, is it ok to treat srid 4326's lat/long points as normal x,y co-ordinates? </div>
<div>I ask because when I try to compute if a point falls within the raster extent or not ... then I can just treat it as trying to figure out if a point lie within a rectangle or not, where I know the lower left and top right corner's (lat,long) co-ordinates. Is that still ok? Or am I overlooking something?</div>
<div><br><br><div class="gmail_quote">On Thu, Mar 22, 2012 at 10:17 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">
The raster was just a quick way to make a vector index grid. If you can't get ST_AsRaster to work by upgrading to GDAL 1.9 use ST_MakeEmptyRaster and then ST_MapalgebraExpr with the '([rast.x] - 1) * ' || ST_Width(rast)::text || ' + [rast.y]' expression.<br>
<br>
Then convert each pixel to a table of polygons with (ST_PixelAsPolygons(rast)).geom & (ST_PixelAsPolygons(rast)).val<br>
<br>
Then proceed with a normal geometry/geometry intersection.<br>
<br>
We don't have a very fast way to do that in raster mode right now. The fastest approach with what we have would probably be to use ST_MapAlgebraFct and write a custom function that generate the values by doing a point/geometry instersection with the road buffer table.<br>
<span class="HOEnZb"><font color="#888888"><br>
Pierre<br>
</font></span><div class="im HOEnZb"><br>
> -----Original Message-----<br>
> From: <a href="mailto:postgis-users-bounces@postgis.refractions.net">postgis-users-bounces@postgis.refractions.net</a> [mailto:<a href="mailto:postgis-users-">postgis-users-</a><br>
> <a href="mailto:bounces@postgis.refractions.net">bounces@postgis.refractions.net</a>] On Behalf Of Ed Linde<br>
</div><div class="HOEnZb"><div class="h5">> Sent: Thursday, March 22, 2012 5:03 PM<br>
> To: PostGIS Users Discussion<br>
> Subject: Re: [postgis-users] ST_Buffer + grid problem<br>
><br>
> Oh and if you think there is a smarter way to do this, like maybe not making a<br>
> raster like I am trying to etc then please let me know. The idea is that the road<br>
> geometry is way way smaller than the point cloud which is in Terabytes, so I am<br>
> making a grid index in some way around the road buffers, to speed up the<br>
> computation.<br>
><br>
><br>
> On Thu, Mar 22, 2012 at 10:01 PM, Ed Linde <<a href="mailto:edolinde@gmail.com">edolinde@gmail.com</a>> wrote:<br>
><br>
><br>
> Hi Pierre,<br>
> The idea is that I want to build this uniform grid on the road geometry<br>
> and like you mention it has dimensions 1020x24798. I want to assign each cell a<br>
> unique cell ID based on some formula like you had earlier ... x*width + y to get a<br>
> 1D cell ID from (x,y). Hopefully I can do this on just (lat.long) of srid 4326<br>
> directly. Then I need to<br>
> intersect it to the buffers around the roads which I think should be easy<br>
> once I have a raster, so this is just a normal intersect. Now I should know each<br>
> cell and what road buffer it belongs to! Then on my second data set which is a<br>
> massive point cloud with srid = 4326 as well, I just want to pass through it ONCE<br>
> and compute in which cell it falls and hence associate a point with a road buffer.<br>
> Hope that made sense :)<br>
><br>
> Cheers,<br>
> Ed<br>
><br>
> On Thu, Mar 22, 2012 at 9:53 PM, Pierre Racine<br>
> <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>> wrote:<br>
><br>
><br>
> 1020x24798 so more than 25 000 000 polygons... or pixels...<br>
><br>
> So if I understand well you want to assign some values to each<br>
> of those cells based on a vector coverage (of how many polygons)?<br>
><br>
><br>
> > -----Original Message-----<br>
> > From: <a href="mailto:postgis-users-bounces@postgis.refractions.net">postgis-users-bounces@postgis.refractions.net</a><br>
> [mailto:<a href="mailto:postgis-users-">postgis-users-</a><br>
> > <a href="mailto:bounces@postgis.refractions.net">bounces@postgis.refractions.net</a>] On Behalf Of Ed Linde<br>
><br>
> > Sent: Thursday, March 22, 2012 4:43 PM<br>
> > To: PostGIS Users Discussion<br>
> > Subject: Re: [postgis-users] ST_Buffer + grid problem<br>
> ><br>
><br>
> > So far no proper raster, because the following failed, but<br>
> maybe the polygon<br>
> > gives us an idea? This polygon is the extent of all the road<br>
> geometries in<br>
> > Denmark. Also I remember that when you did the raster at<br>
> your end you said it<br>
> > worked. Is 0.000036 degrees = 4m correct if I want to get a<br>
> 4m by 4m cell sized<br>
> > uniform raster grid?<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>
> ><br>
> > On Thu, Mar 22, 2012 at 9:39 PM, Pierre Racine<br>
> <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>><br>
> > wrote:<br>
> ><br>
> ><br>
> > > Thanks. The raster I need to visualise is a 4m by 4m grid<br>
> on the entire<br>
> > map of<br>
> > > Denmark! :) So do you classify that as a large raster? If so<br>
> is there a<br>
> > way to see a<br>
> > > portion of it or something?<br>
> > > I just want to manually check for a few areas on the map<br>
> that the<br>
> > intersection<br>
> > > indeed works as expected and things haven't gone awry<br>
> thanks to<br>
> > projection<br>
> > > differences or that I had degrees instead of meters or<br>
> some such<br>
> > thing. Well, I<br>
> > > still have to first test if this GDAL upgrade will fix things<br>
> and if make<br>
> > empty raster<br>
> > > .. makes a difference.<br>
> ><br>
> ><br>
> > That must be big but Danemark is a small country ;-) Only<br>
> width and<br>
> > height tell you if a raster is big.<br>
> ><br>
> > In the database, your grid is a set of rectangular geometry<br>
> (how many?)<br>
> > or just a big raster now (width & height)?<br>
> ><br>
> > _______________________________________________<br>
> > postgis-users mailing list<br>
> > <a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br>
> > <a href="http://postgis.refractions.net/mailman/listinfo/postgis-" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-</a><br>
> users<br>
> ><br>
> ><br>
><br>
> _______________________________________________<br>
> postgis-users mailing list<br>
> <a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><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>
postgis-users mailing list<br>
<a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br>
<a href="http://postgis.refractions.net/mailman/listinfo/postgis-users" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a><br>
</div></div></blockquote></div><br></div>