rasterize_lidar_cover_class() ignores low vegetation

[delphinekw]

Hi and thanks for the wonderful package!

I have been trying to create shade maps for a Swiss city. The lidar data that is published by the Swiss government uses only one class for vegetation: low vegetation.

The function rasterize_lidar_cover_class() has an inbuilt correction of low vegetation through following line, which reclassifies low vegetation to ground points:
pts$Classification[pts$Classification == 3] <- as.integer(2)

In my case, this creates incomplete rasters for cover classes, and the resulting shade raster at ground level is incorrect.

I would suggest removing this line from the function definition, or alternatively, making it optional with an extra argument in the function definition, i.e. correct.low.veg = T that enables deactivating this behaviour.

Thanks again for a very useful package!

[delphinekw]

I just created a pull request with the suggested changes. I hope this helps

[Pakillo]

Hi @delphinekw

Thanks for reporting and the suggestions. I could add the logical argument to apply the correction, but I'm unsure how that would affect the rest of the package, which assumes that pixels are either (i) buildings, (ii) high vegetation (which can cast shade on the ground) and (iii) ground, including low vegetation that doesn't cast shade at a practical level (that's why low vegetation is merged with ground).

The function rasterize_lidar_cover_class() is there to be able to close the entire workflow with CityShadeMapper, but it can be fragile with some lidar point clouds. I guess you tried with fill.holes = TRUE, right? That manages to solve some issues. You could also do some post-processing of the resulting raster to fix issues. But if you can obtain a raster of ground, tall vegetation, and buildings cover by some other means (e.g. from the city council or high-resolution land cover raster) that may be more reliable. You could also try producing such a raster from image segmentation (e.g. see https://github.com/walkerke/geosam).

Anyway, I've opened a branch with a new argument low.veg.as.ground following your suggestion. You can install it using
remotes::install_github("Pakillo/CityShadeMapper", ref = remotes::github_pull("23")). Please let me know how it works.

I'm preparing a funding application to support further development of this package and I'd be very grateful if you could drop me an email (see https://frodriguezsanchez.net/#contact) about your use case and other needs or features you'd like to see implemented. I'd try to get them in further development plans. I want this software to be useful and help cities' environmental planning. Thanks

[Pakillo]

The lidar data that is published by the Swiss government uses only one class for vegetation: low vegetation.

Ah, I'd missed this. I understand all vegetation is then classified as low in that lidar dataset, regardless their height. The fix in #23 should work fine then, as vegetation height is considered apart in the process. But if you want to calculate shade at the ground (rather than canopy level), note that the ground illumination level under vegetation is going to be fixed at 5%. Which may make sense for trees, but maybe not so much for low vegetation < 1m high. A person walking over that terrain is going to receive much more illumination than 5% unless there are taller trees or buildings casting some shade. I'd suggest to reclassify the ground illumination values conditional on vegetation height afterwards. HTH

[delphinekw]

Thanks for your quick reply!

Yes exactly, all the vegetation is classified under low vegetation in the swiss public datasets (at least in my area). And there are no other available datasets for smaller municipalities.

I see your point. Your suggestion to reclassify the ground illumination based on vegetation height makes sense, I'm just a little unsure how that would work in practice. There might be some induced errors because the angle of the sunlight is taken into consideration in the ground illumination but not the vegetation height raster, especially if you have low vegetation right next to high vegetation.

Probably it would then be better to reclassify my point cloud to separate low and high vegetation. For this, I would need to dive back into the lidR package. Let's see when I find the time :-)

Another fast workaround I tried was to reclassify the point cloud when loading for example with a filter from the lidR package: opt_filter(ctg) <- "-change_classification_from_to 3 4", but this would need to be permitted in the rasterize_lidar_cover_class() function, that itself already defines some filters (it might anyway be a good idea to enable the use of other filters, e.g. -keep_xy to test the functions on smaller areas)

[delphinekw]

The changes in the new branch work well for me. I obtain ground illumination maps that make sense.

The quality issue you mentioned remains (the shade cast by low vegetation is useless for a pedestrian). But for a first version of the shade maps, and to brainstorm how we can use them, it will work well.

I will get in touch, as soon as I have discussed some applications with my team. Thank you very much for your time!

[Pakillo]

Probably it would then be better to reclassify my point cloud to separate low and high vegetation.

Yes, agree this looks like the best approach.

I think you could easily redefine all vegetation points above a certain height (e.g. 2 meters) like this:
las$Classification[las$Classification == 3 & las$Z > 2] <- as.integer(4)

Full reproducible example:

library(CityShadeMapper)

# read las
example <- system.file("extdata", "PlazaNueva.laz", package = "CityShadeMapper")
las <- lidR::readLAS(example, select = "xyzc")

table(las$Classification)
#> 
#>     1     2     3     4     5     6     7    12 
#>    81 11935  2784   697 17137 18842    72 51017

# reclassify all low vegetation points with height >2m as medium vegetation (class 4)
las$Classification[las$Classification == 3 & las$Z > 2] <- as.integer(4)

# result
table(las$Classification)
#> 
#>     1     2     4     5     6     7    12 
#>    81 11935  3481 17137 18842    72 51017

^{Created on 2026-01-28 with reprex v2.1.1}

Then you could safely mix low vegetation together with ground cover, and keep vegetation >2m tall to account for its shading.

[delphinekw]

Thank you very much for your example!

I am still struggling a bit with the data handling between CityShadeMapper and lidR. Your example will change the data in RAM but not on the drive, and the function rasterize_lidar_cover_class() will load the data from the drive again, and the changes will be lost.

I can work around this, and loop over my LAS tiles and save them to the drive, but I feel like there could be a better solution. What am I missing?

What would be the "proper" way of using CityShadeMapper with your reproducible example? The last line in this example throws an error:

`library(CityShadeMapper)

example <- system.file("extdata", "PlazaNueva.laz", package = "CityShadeMapper")
las <- lidR::readLAS(example, select = "xyzcr")
las$Classification[las$Classification == 3 & las$Z > 2] <- as.integer(4)

heights <- calc_heights_from_lidar(las)

cover.ras <- rasterize_lidar_cover_class(las)
#> Fehler in UseMethod("readLAS", files) :
nicht anwendbare Methode für 'readLAS' auf Objekt der Klasse "LAS" angewendet`

[Pakillo]

Hi

I've modified the function to allow taking LAS objects existing in memory, rather than always reading them from disk. So your example should work now. Please reinstall the package from github and let me know

Cheers

[delphinekw]

Hi

It throws a different error: "Object 'pts' not found"

Thanks

[Pakillo]

sorry, a typo, should be fixed now

[delphinekw]

Works like a charm! Thank you for your support

Maybe this makes the modifications you made yesterday in #23 less relevant, since there is a good reason to ignore low vegetation. But I will leave that jugdment to you :-)