Thanks everyone for the help and tips. I ended up using a tool called Lidar Fusion, from the forestry people at USDA, USFS and university of washington. http://forsys.cfr.washington.edu/
It works pretty fast and does exactly what I needed. But only found executables for windows... Cheers Daniel On Fri, Mar 30, 2012 at 12:22 PM, Daniel Victoria <[email protected]> wrote: > Werner, > > Laslib and specially lasheight looks like what I want. However, I > already have the ground raster interpolated by the people that > provided the las data so I'll have to check if I can ignore their > ground raster and use the one generated by laslib. But many thanks for > the help, I'll certainly study that option more. > > Michael, > > I liked your python suggestion. Will give it a try. I'll just have to > study your script a little more in order to figure out if it's > considering the case where there are more than 1 point per ground > cell. > > My ground raster has 1m resolution but my point density is at least 4 > points per meter... > > Cheers to all > Daniel > > On Fri, Mar 30, 2012 at 11:52 AM, Werner Macho <[email protected]> wrote: >> Hi! >> Don't know if i understood everything correctly. (Only ob mobile right now) >> If it's not mandatory to use grass you can also use Laslib drom martin >> isenburg. I am pretty sure that. it should do what you want pretty fast and >> with huge amount of points.. >> I hope i got you right .. if not sorry f0r the noise >> >> Regards >> Werner >> >> Am 30.03.2012 16:43 schrieb "michael vetter" <[email protected]>: >> >>> hey, >>> >>> you can also calculate a normalized point cloud by subtract Z(from point) >>> from the DTM(Z) >>> >>> you need a DTM (as raster) and the LiDAR points as ASCII >>> >>> >>> 1. "r.stats -gn in=DTM >temp_elev" >>> 2. create a dictionary in python, from "temp_elev" k=x,y and value=z (x >>> and y should INT) >>> 3. loop the point cloud and subtract the z from k for each LiDAR point and >>> add a nZ column >> nDTM_LiDAR_points >>> >>> 2 and 3 should run in the same script! >>> >>> e.g. In Python: >>> >>> #read input raster data (ascii from r.stats -gn DTM) >>> for lineR in iR.readlines(): >>> valR = lineR.strip().split(' ') >>> x = float(valR[0]) >>> y = float(valR[1]) >>> k = ('%s %s' % (int(x),int(y))) >>> z = float(valR[2]) >>> value.append(z) >>> D[k] = value >>> value = [] >>> >>> #read input point cloud as ascii and caclate nZ >>> for lineP in iP.readlines(): >>> valP = lineP.strip().split(' ') >>> x = float(valP[0]) >>> y = float(valP[1]) >>> k = ('%s %s' % (int(x),int(y))) >>> Z = float(valP[2]) >>> output.write('%s %s %s %s\n' %(x,y,Z,Z-D[k][0])) >>> oF.write(output.getvalue()) >>> >>> then you can 'awk' the data you are interested in 'cat nDTM_LiDAR_points | >>> awk '{if($4 < 2 && $4 > 1) print $0}' >LiDAR_1_to_2 >>> >>> r.in.xyz LiDAR_1_to_2 out=LiDAR_1_to_2 method=n >>> >>> you can see the related reference to a similar workflow: >>> >>> http://www.isprs.org/proceedings/XXXVIII/5-W12/Papers/ls2011_submission_35.pdf >>> >>> Michael >>> >>> Am 2012-03-30 14:52, schrieb Daniel Lee: >>> >>> Hi Daniel V ;) >>> >>> You're right, that might not be the best way to go. I thought that it >>> might simply be faster to do a topological operation rather than a DB edit. >>> To be honest, I'd stay away from any 3D objects like volumes because it'd >>> just get pretty complicated if you use them... As far as I know ;) Could be >>> wrong. My suggestion would use the following steps: >>> >>> 1. Make terrain raster >>> 2. Make surface raster >>> -- Here I'm assuming that for you basically have only two height "layers" >>> - i.e. no points that contain information between the surface raster and >>> terrain (like bushes beneath a tree cover). >>> 3. Make a relative digital surface model (rDSM) (r.mapcalc --> rDSM = >>> surface - terrain) >>> 4. Reclassify the rDSM into the different classes you're interested in >>> 5. Convert the reclassified raster into an area vector file. >>> -- This makes a 2D map of polygons that cover the same areas as the raster >>> classes. >>> 6. Extract the points inside the different different height classes >>> (v.select) >>> >>> However, if you're working with regularly spaced points that pass onto a >>> regular grid (=1 pt. / raster pixel) you could always use r.stats to tell >>> you how much area is in each raster category. Then you really wouldn't need >>> vectors at all. >>> >>> Daniel L. >>> >>> -- >>> >>> B.Sc. Daniel Lee >>> Geschäftsführung für Forschung und Entwicklung >>> ISIS - International Solar Information Solutions GbR >>> Vertreten durch: Daniel Lee, Nepomuk Reinhard und Nils Räder >>> >>> Deutschhausstr. 10 >>> 35037 Marburg >>> Festnetz: +49 6421 379 6256 >>> Mobil: +49 176 6127 7269 >>> E-Mail: [email protected] >>> Web: http://www.isi-solutions.org >>> >>> ISIS wird gefördert durch die Bundesrepublik Deutschland, Zuwendungsgeber: >>> Bundesministerium für Wirtschaft und Technologie aufgrund eines Beschlusses >>> des Deutschen Bundestages, sowie durch die Europäische Union, >>> Zuwendungsgeber: Europäischer Sozialfonds. >>> Zusätzliche Unterstützung erhält ISIS von dem Entrepreneurship Cluster >>> Mittelhessen, der Universität Marburg, dem Laboratory for Climatology and >>> Remote Sensing und dem GIS-Lab Marburg. >>> >>> >>> >>> >>> Am 30. März 2012 14:09 schrieb Daniel Victoria >>> <[email protected]>: >>>> >>>> Hi Daniel Lee (two Daniels changing emails is a hard thing to follow...) >>>> >>>> I can create the different height classes raster easily but I'm not >>>> sure how to get the lidar points from the cloud that are in each >>>> vertical slice. >>>> >>>> Maybe work with raster volumes? Is there a way to cross points at >>>> different volumes? >>>> >>>> The fact is, I'm not seeing how to do this without having to import >>>> the lidar point cloud. >>>> >>>> Schematically, what I want to do is count the number of points that >>>> are in the vertical bins 1 2, (g is ground), a vertical slice... >>>> ___ >>>> | . : >>>> bin 2 | : >>>> | . >>>> ----- >>>> bin 1 | ... >>>> | . >>>> ggggg >>>> >>>> Cheers and many thanks for the attention >>>> Daniel V >>>> >>>> >>>> On Fri, Mar 30, 2012 at 8:08 AM, Daniel Lee <[email protected]> >>>> wrote: >>>> > Hi Daniel, >>>> > >>>> > You could try making different rasters with classes (0-5m over ground, >>>> > 5-10m >>>> > over ground, etc.) and then convert them into polygons, then check how >>>> > many >>>> > points are inside them. I'm not sure if it'd be faster than the way you >>>> > suggested originally, though, when you think that you'd have to make >>>> > the >>>> > rasters first, etc. >>>> > >>>> > I ended up giving up on large vector operations with LiDAR point clouds >>>> > in >>>> > GRASS 6.4 because it simply took way too long, but I was dealing with >>>> > millions of points. For such a small dataset it's probably fine, since >>>> > as >>>> > long as you know your script is okay you can let it run through your >>>> > lunch >>>> > break ;) >>>> > >>>> > >>>> > Daniel >>>> > >>>> > -- >>>> > >>>> > B.Sc. Daniel Lee >>>> > Geschäftsführung für Forschung und Entwicklung >>>> > ISIS - International Solar Information Solutions GbR >>>> > Vertreten durch: Daniel Lee, Nepomuk Reinhard und Nils Räder >>>> > >>>> > Deutschhausstr. 10 >>>> > 35037 Marburg >>>> > Festnetz: +49 6421 379 6256 >>>> > Mobil: +49 176 6127 7269 >>>> > E-Mail: [email protected] >>>> > Web: http://www.isi-solutions.org >>>> > >>>> > ISIS wird gefördert durch die Bundesrepublik Deutschland, >>>> > Zuwendungsgeber: >>>> > Bundesministerium für Wirtschaft und Technologie aufgrund eines >>>> > Beschlusses >>>> > des Deutschen Bundestages, sowie durch die Europäische Union, >>>> > Zuwendungsgeber: Europäischer Sozialfonds. >>>> > Zusätzliche Unterstützung erhält ISIS von dem Entrepreneurship Cluster >>>> > Mittelhessen, der Universität Marburg, dem Laboratory for Climatology >>>> > and >>>> > Remote Sensing und dem GIS-Lab Marburg. >>>> > >>>> > >>>> > >>>> > >>>> > Am 30. März 2012 13:00 schrieb Daniel Victoria >>>> > <[email protected]>: >>>> > >>>> >> I need to calculate the number of points at different height levels. >>>> >> That >>>> >> is, how many points are from 0 to 5 meters? And from 5 to 10? And so >>>> >> on. So >>>> >> I believe I need to work with vector points and a database. Or is >>>> >> there >>>> >> another way? >>>> >> Thanks >>>> >> Daniel >>>> >> >>>> >> On Mar 30, 2012 6:47 AM, "Daniel Lee" <[email protected]> wrote: >>>> >>> >>>> >>> Hi there, >>>> >>> >>>> >>> I don't work too much with vectors, but my personal experience with >>>> >>> them >>>> >>> has been fairly slow, perhaps due to the topology. If you have access >>>> >>> to the >>>> >>> raw point cloud, you could try importing the ground and surface >>>> >>> points as >>>> >>> separate rasters using r.in.xyz and then use r.mapcalc to get the >>>> >>> height by >>>> >>> subtracting the ground from the surface raster. Or do you definitely >>>> >>> need >>>> >>> vector points as an output? >>>> >>> >>>> >>> Best, >>>> >>> Daniel >>>> >>> >>>> >>> -- >>>> >>> >>>> >>> B.Sc. Daniel Lee >>>> >>> Geschäftsführung für Forschung und Entwicklung >>>> >>> ISIS - International Solar Information Solutions GbR >>>> >>> Vertreten durch: Daniel Lee, Nepomuk Reinhard und Nils Räder >>>> >>> >>>> >>> Deutschhausstr. 10 >>>> >>> 35037 Marburg >>>> >>> Festnetz: +49 6421 379 6256 >>>> >>> Mobil: +49 176 6127 7269 >>>> >>> E-Mail: [email protected] >>>> >>> Web: http://www.isi-solutions.org >>>> >>> >>>> >>> ISIS wird gefördert durch die Bundesrepublik Deutschland, >>>> >>> Zuwendungsgeber: Bundesministerium für Wirtschaft und Technologie >>>> >>> aufgrund >>>> >>> eines Beschlusses des Deutschen Bundestages, sowie durch die >>>> >>> Europäische >>>> >>> Union, Zuwendungsgeber: Europäischer Sozialfonds. >>>> >>> Zusätzliche Unterstützung erhält ISIS von dem Entrepreneurship >>>> >>> Cluster >>>> >>> Mittelhessen, der Universität Marburg, dem Laboratory for Climatology >>>> >>> and >>>> >>> Remote Sensing und dem GIS-Lab Marburg. >>>> >>> >>>> >>> >>>> >>> >>>> >>> >>>> >>> Am 29. März 2012 22:59 schrieb Daniel Victoria >>>> >>> <[email protected]>: >>>> >>>> >>>> >>>> Hi all, >>>> >>>> >>>> >>>> I'm trying to calculate the height from the ground of several lidar >>>> >>>> points (15million) in order to get the number of points that occur >>>> >>>> at >>>> >>>> different height levels. I read some older posts about using >>>> >>>> r.in.xyz >>>> >>>> (or r.in.lidar in grass 7) but I could not understand how to count >>>> >>>> the >>>> >>>> number or points that have height from 0 to 5 m above the ground, >>>> >>>> for >>>> >>>> instance. So, I'm trying to do the following. >>>> >>>> >>>> >>>> 1) import points using v.in.lidar (Grass 7 ubuntu linux) >>>> >>>> 2) create a column in the database for the ground height and one for >>>> >>>> elevation >>>> >>>> 3) populate height column from ground raster (generate in another >>>> >>>> process) using v.what.rast >>>> >>>> 4) calculate elevation as zcoord - ground for each point >>>> >>>> (v.db.update) >>>> >>>> >>>> >>>> As you might imagine, this takes a long time. Just to import a 400Mb >>>> >>>> lidar file takes around 50min. >>>> >>>> Is there any easier way that I'm not envisioning? >>>> >>>> >>>> >>>> I'm running Grass 7 with liblas on a Ubuntu Virtual Machine and >>>> >>>> sqlite >>>> >>>> backend. >>>> >>>> >>>> >>>> Thanks >>>> >>>> Daniel >>>> >>>> _______________________________________________ >>>> >>>> grass-user mailing list >>>> >>>> [email protected] >>>> >>>> http://lists.osgeo.org/mailman/listinfo/grass-user >>>> >>> >>>> >>> >>>> > >>> >>> >>> >>> >>> _______________________________________________ >>> grass-user mailing list >>> [email protected] >>> http://lists.osgeo.org/mailman/listinfo/grass-user >>> >>> >>> >>> -- >>> Mag. Michael Vetter >>> >>> Institute of Photogrammetry and Remote Sensing (I.P.F.) >>> Vienna University of Technology (TU Wien) >>> Gusshausstrasse 27-29, 1040 Vienna, Austria >>> >>> >>> Centre for Water Resource Systems (CWRS) >>> Vienna University of Technology (TU Wien) >>> Karlsplatz 13/222, A-1040 Vienna, Austria >>> >>> >>> Tel: +43-(0)1-58801-22226 >>> E-mail: [email protected] >>> www.ipf.tuwien.ac.at >>> www.waterresources.at >>> >>> >>> >>> >>> _______________________________________________ >>> grass-user mailing list >>> [email protected] >>> http://lists.osgeo.org/mailman/listinfo/grass-user >>> >> >> _______________________________________________ >> grass-user mailing list >> [email protected] >> http://lists.osgeo.org/mailman/listinfo/grass-user >> _______________________________________________ grass-user mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-user
