I need to obtain a grid of bare ground elevation (removing all edifices, trees, roads and man-made elements) for quantitative geomorphological analysis from a high resolution (1 m) DSM obtained by orthoimages (not from LIDAR). Is there any suggestion you can provide me to filter my original grid?
DEM from DSM – How to Obtain a DTM/DEM from DSM
demelevation
Related Solutions
In my experience, DEM is most of the time used as a generic term for DSMs and DTMs. I think this image on Wikipedia depicts the differences between DSMs and DTMs well:
- DSM = (earth) surface including objects on it
- DTM = (earth) surface without any objects
A different definition is found in [Li et al., DIGITAL TERRAIN MODELING - Principles and Methodology]:
DEM is a subset of DTM and the most fundamental component of DTM.
In practice, these terms (DTM, DEM, DHM, and DTEM) are often assumed to be synonymous and indeed this is often the case. But sometimes they actually refer to different products. That is, there may be slight differences between these terms. Li (1990) has made a comparative analysis of these differences as follows:
- Ground: “the solid surface of the earth”; “a solid base or foundation”; “a surface of the earth”; “bottom of the sea”; etc.
- Height: “measurement from base to top”; “elevation above the ground or recognized level, especially that of the sea”; “distance upwards”; etc.
- Elevation: “height above a given level, especially that of sea”; “height above the horizon”; etc.
- Terrain: “tract of country considered with regarded to its natural features, etc.”; “an extent of ground, region, territory”; etc.
R and Fortran will read an ArcInfo binary grid directly, with essentially no extra processing, so skip the middleman: you don't need QGIS or anything else.
These datasets consist of two files. One is an ASCII header file with .hdr extension, formatted to look like this:
ncols 1133
nrows 1415
xllcorner 686280
yllcorner 4179990
cellsize 30
NODATA_value -9999
byteorder LSBFIRST
The meanings should be evident.
The other file is an array of four-byte values (either floats or integers, depending on the type of the grid) arranged row by row, from top down.
Here is a quick and dirty example of how it can be read and processed in R. (It's quick and dirty because it does not check for deviations from the format or inconsistencies, has no error checking, and has not been tested on all possible systems.) To demonstrate that the data manipulation is correct, I plot the data and show their coordinates.
Automation is trivial, of course: encapsulate this within a function and enclose it within an appropriate iterator (*apply, for, while, or whatever in R; FOR, DO, WHILE, or whatever in Fortran).
R's raster package supports subsetting of grids by general polygons. For rectangular subsetting, that's a basic part of R's array indexing capabilities.
#
# Read the header file.
#
sFile <- "G:/USGS/DEM/7_5min/VA/albem_s" # Root dataset name, without extensions
h <- read.table(paste(sFile, ".hdr", sep=""),
col.names=c("field", "value"), as.is=TRUE)
#
# Extract the header values.
#
header <- type.convert(h[1:6,"value"]); names(header) <- h[1:6,"field"]
n.rows <- header["nrows"]
n.cols <- header["ncols"]
nodata <- header["NODATA_value"]
x0 <- header["xllcorner"]
y0 <- header["yllcorner"]
cellsize <- header["cellsize"]
byteorder <- ifelse(h[7, "value"]=="LSBFIRST", "little", "big")
#
# Read the data.
#
x <- readBin(paste(sFile, ".flt", sep=""),
"numeric", n=n.rows*n.cols, size=4, endian=byteorder)
x[x==nodata] <- NA # Mark the NoData values
x <- t(matrix(x, nrow=1133)) # (R stores data by columns, not rows)
#
# Display the grid.
#
library(raster)
plot(raster(x, xmn=x0, xmx=x0+n.cols*cellsize, ymn=y0, ymx=y0+n.rows*cellsize))
This image looks correct. The total time to read and process the dataset (before displaying it) was 1/8 second.
Best Answer
This filtering process is usually performed on a lidar point cloud and not an interpolated derivative. It is unlikely that you will have satisfactory results attempting to filter the DSM. I would highly recommend tracking down the original lidar data.
You could attempt to treat your DSM as a point cloud by converting it to points and then running a filter intended for lidar point clouds. Depending on the algorithm, you may get a suitable result for generating a bare earth DEM. However, it may end up oversmoothed and not supporting the current resolution of your DSM.
Some recommended "free" lidar filtering software:
Airborne LIDAR Data Processing and Analysis Tools (ALDPAT)
GRASS GIS specifically, v.lidar.correction
Idaho State University Boise Center Aerospace Lab IDL Virtual Machine software (BCAL)
SAGA GIS
USFS-PNW lidar processing and visualization software (FUSION)
USFS-RMRS Multiscale Curvature Classification (MCC)