[GRASS-user] Lidar points height from ground
Daniel Victoria
daniel.victoria at gmail.com
Fri Mar 30 16:33:59 EDT 2012
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
>>
More information about the grass-user
mailing list