[GRASS-user] Lidar points height from ground
michael vetter
mv at ipf.tuwien.ac.at
Fri Mar 30 10:28:01 EDT 2012
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 <mailto:Lee at isi-solutions.org>
> Web: http://www.isi-solutions.org <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 <mailto: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
> <mailto: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 <tel:%2B49%206421%20379%206256>
> > Mobil: +49 176 6127 7269 <tel:%2B49%20176%206127%207269>
> > E-Mail: Lee at isi-solutions.org <mailto: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 <mailto: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
> <mailto: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 <tel:%2B49%206421%20379%206256>
> >>> Mobil: +49 176 6127 7269 <tel:%2B49%20176%206127%207269>
> >>> E-Mail: Lee at isi-solutions.org <mailto: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 <mailto: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 <mailto: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
-------------- nächster Teil --------------
Ein Dateianhang mit HTML-Daten wurde abgetrennt...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20120330/24927d35/attachment-0001.html
More information about the grass-user
mailing list