[GRASS-user] Lidar points height from ground

Sat Mar 31 07:03:08 EDT 2012

Hi there,

Bit late in the game to this discussion but the tool I use for all my lidar point processing (on Windows and linux) is http://www.cs.unc.edu/~isenburg/lastools/

I've found it's very powerful and user friendly (but still "tunable" to your preferences) and the users list is also really good...

Happy processing!

Rebecca

Dr Rebecca Bennett
Researcher (Archaeology Group)
School of Applied Sciences, Bournemouth University

>________________________________
> From: Daniel Victoria <daniel.victoria at gmail.com>
>To: Werner Macho <werner.macho at gmail.com> 
>Cc: grass-user at lists.osgeo.org 
>Sent: Friday, 30 March 2012, 21:33
>Subject: Re: [GRASS-user] Lidar points height from ground
> 
>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
><daniel.victoria at gmail.com> 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 <werner.macho at gmail.com> 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" <mv at ipf.tuwien.ac.at>:
>>>
>>>> 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: Lee at isi-solutions.org
>>>> 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
>>>> <daniel.victoria at gmail.com>:
>>>>>
>>>>> 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 <lee at isi-solutions.org>
>>>>> 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: Lee at isi-solutions.org
>>>>> > 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
>>>>> > <daniel.victoria at gmail.com>:
>>>>> >
>>>>> >> 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" <lee at isi-solutions.org> 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: Lee at isi-solutions.org
>>>>> >>> 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
>>>>> >>> <daniel.victoria at gmail.com>:
>>>>> >>>>
>>>>> >>>> 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
>>>>> >>>> grass-user at lists.osgeo.org
>>>>> >>>> http://lists.osgeo.org/mailman/listinfo/grass-user
>>>>> >>>
>>>>> >>>
>>>>> >
>>>>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> grass-user mailing list
>>>> grass-user at lists.osgeo.org
>>>> 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: mv at ipf.tuwien.ac.at
>>>> www.ipf.tuwien.ac.at
>>>> www.waterresources.at
>>>>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> grass-user mailing list
>>>> grass-user at lists.osgeo.org
>>>> http://lists.osgeo.org/mailman/listinfo/grass-user
>>>>
>>>
>>> _______________________________________________
>>> grass-user mailing list
>>> grass-user at lists.osgeo.org
>>> http://lists.osgeo.org/mailman/listinfo/grass-user
>>>
>_______________________________________________
>grass-user mailing list
>grass-user at lists.osgeo.org
>http://lists.osgeo.org/mailman/listinfo/grass-user
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20120331/0b903ce5/attachment-0001.html