[GRASS-user] Lidar points height from ground

Daniel Victoria daniel.victoria at gmail.com
Mon Apr 2 08:22:25 EDT 2012 

Hi al,

Very good to know about lastools. Will certainly give it a try.

Cheers
Daniel

On Sun, Apr 1, 2012 at 4:31 PM, Martin Isenburg
<martin.isenburg at gmail.com> wrote:
> Hello,
>
> thanks Rebecca for the kind words. With the most recent release of
> LAStools you can also use lasheight to compute points above an
> existing raster DTM givem that this raster DTM is in ASC or XYZ
> format. A similar question was just asked on the LAStools user forum.
> For the full thread see here:
>
> http://groups.google.com/group/lastools/browse_thread/thread/6159fbeb1b9af720
>
> The answer was to make use of the fact that as of the 120327 release
> all LAStools can read ESRI's ASCII rasters *.asc as if they were LiDAR
> point files. So simply convert your current raster DEM to the *.asc
> format and run
>
> lasheight -i lidar.laz -ground_points dem.asc -replace_z -o lidar_above_dem.laz
>
> cheers,
>
> Martin @lastools
>
> --
> http://lastools.org - http://laszip.org
>
> On Sat, Mar 31, 2012 at 4:03 AM, Rebecca Bennett <rabennett at ymail.com> wrote:
>> 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
>>
>>
> _______________________________________________
> grass-user mailing list
> grass-user at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user

More information about the grass-user mailing list