[GIS] Obtain all elevations points from a raster DEM

demelevationqgis

I need to obtain all elevation grid points from a raster DEM.

I see the steps required as:

  1. Read in raster DEM (Arc Info Binary Grid format).
  2. Create subset area using Points2One (see triangle area here: http://imgur.com/aeOPBBD).
  3. Write subset-grid elevations to a file (I can't do this yet). I only need elevations (ie. Z only, not XY).

Any ideas for step 3? I'm happy with getting all grid elevations or a histogram of elevations (with frequencies for (say) 100 bins). I will later manipulate the elevation data in EXCEL or R or Fortran or Python or ???

Note: I'm very new to using GISs and am using QGIS 1.8.0 on Windows (since I found that the easiest program to use). Thus, I'd prefer a solution using QGIS, but I'm willing to use another free GIS program if required. Also, I'll likely have to automate/program this extraction procedure later (ie. steps 2 and 3), but would prefer a GUI-based solution first for simplicity.

Thanks, Tom.

Best Answer

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))

Figure

This image looks correct. The total time to read and process the dataset (before displaying it) was 1/8 second.

Related Question