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
