[GRASS-user] Lidar points height from ground
Daniel Victoria
daniel.victoria at gmail.com
Fri Mar 30 11:22:25 EDT 2012
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