I would like to use 3D spatial metrics that introduced in the paper titled "Lidar and VHRS data for assessing living quality in cities - An approach based on 3D spatial indices". I have a lidar dataset including 14 las files and I have a shapefile including 416 buildings. I wanted to calculate cubic volume of each building using LiDAR data. As I understand,I need to clip las catalog for each object in the shapefile and then run voxels on every clipped las file. I tried to do that but I could not understand how to get the volume calculation for each building. Is there any way to write volume value into the shapefile?

  • Please edit your question to explain what shp contains. A print at least + maybe an image. And clean you question to remove code that is not relevant. I see nothing that looks like an attempt to measure building volume and I see many lines that are not related at all to the question. Also explain how you want to compute building volume. Is the point cloud classified? Is your shapefile the reference?
    – JRR
    Mar 26, 2020 at 10:45
  • Thanks, I edited my question. Unfortunately, I do not have any print except my codes yet Mar 26, 2020 at 18:14

1 Answer 1


Assuming that your shapefile contains polygons you could loop through each polygon and extract the points within the polygons.

For each point cloud, it is expected to contain only one building if your shapefile is well build. Also you have the area of the building with the polygon.

Thus, you can analyses each point cloud to find the elevation of your building. If we assume that the buildings are cuboid the highest point can do the job. Here what it could look like:

ctg = readLAScatalog(...)
buildings = shapefile(...)
buildings$Volume = 0

for (i in 1:length(buildings) {
  building = buildings[i,] # Get the polygon i
  las = clip_roi(ctg, building)
  hmax = max(las$Z) # Assuming a normalized dataset
  A = rgeos::gArea(building)
  V = A*hmax
  building$Volume[i] = V

shapefile(buildings, "/where/to/write/building.shp")

If building are not cuboid and you want something more accurate try to get only the surface points and average their heights

las = clip_roi(ctg, building)
las = filter_surfacepoints(las, 2)
hmean = mean(las$Z) # Assuming a normalized dataset
A = rgeos::gArea(building)
V = A*hmean

If you want something yet more accurate you have enough material to tweak these examples.

  • 1
    Thanks a lot. I will never forget this big support! Mar 26, 2020 at 18:38
  • oh I got an Error: unexpected '{' in "for (i in 1:length(buildings) {" > buildings = buildings[i,] # Get the polygon i Error in buildings[i, ] : object 'i' not found > las = lasclip(ctg, buildings) Mar 26, 2020 at 19:21
  • I wrote this code in the forum form. It was not tested. It it just to help you.
    – JRR
    Mar 26, 2020 at 19:23
  • I know, I change as building just in case. Because both of them gave the same error. Mar 26, 2020 at 19:26
  • for (i in 1:length(buildings) { Error: unexpected '{' in "for (i in 1:length(buildings) {" > buiding = buildings[i,] # Get the polygon i Error in buildings[i, ] : object 'i' not found > las = lasclip(ctg, building) Error in lasclip(ctg, building) : object 'building' not found > hmax = max(las$Z) # Assuming a normalized dataset Error: object 'las' not found > A = rgeos::gArea(building) Error in is.projected(spgeom) : object 'building' not found > V = A*hmax Error: object 'A' not found > building$Volume[i] = V Error: object 'V' not found > } Mar 26, 2020 at 19:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.