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.

  • 2
    Because your DEM is already in ArcInfo Binary Grid format, it can immediately be read using Fortran, so there's nothing left to do. (It's just a raw list of four-byte numbers laid out in successive rows; a single Fortran READ statement will suck the whole thing into a 2D array.) – whuber May 21 '13 at 15:33

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

  • 3
    I think you scared him. – nickves May 21 '13 at 16:38
  • Maybe a little scared, but that R seems reasonably simple to me. Unfortunately, my "hdr.adf" seems to be binary, not ASCII - the file starts with "GRID1.2" then some non-ASCII chars (the docs say: "Stored data format: DIGITAL - ArcGIS-grid ArcInfo grid") . Anyway, for now Gene's XYZ solution seems to work for me - although I may investigate an R-based solution later on. Thanks. – tom m May 22 '13 at 3:51
  • With computers, Tom, you have to read things literally: a file named "hdr.adf" is not a .hdr file! You are looking at an ArcInfo grid, rather than a binary export. R will read that directly (with one line) using, for example, the gdal add-in. – whuber May 22 '13 at 14:06

And also with GRASS GIS, SAGA GIS, alone or with the Sextante plugin in QGIS, or R with with other geospatial packages, even with Python (with GDAL/OGR module), or GDAL from the command line, see Clipping rasters with GDAL using polygons but also only with QGIS:

To extract the z values after 2) Go to Raster menu,Conversion sub menu,convert, and choose, for example: (you have the list of all the possible formats in GDAL raster Format)

  • [GDAL]XYZ -- ASCII Gridded XYZ (X,Y,Z format): the result is a text file with the x,y coordinates of the center of the cells of the grid,and the value of the cells (z). After that, you can eliminate the x,y coordinates in the resulting file.

  • [GDAL]ArcInfo -- ASCII Grid(.asc, .ASC) gives you only the z values in a grid form (Esri grid, line,column). If you can program, Python for example, you can extract these z values.

For cropping the image from the command line with GDAL, see Clipping rasters with GDAL using polygons, for example.

  • Thanks for the mention, Gene, but I want to clarify that my solution reads the binary grid directly. It does not propose converting from binary to ASCII and then reading the ASCII grid. It seems you have merely replaced one problem with another that is of equal difficulty (or ease): either read the binary format or read the ASCII format. It involves the same coding either way (although it's arguably now a little easier to access in Excel, fwiw). So the conversion to ASCII just seems like an unnecessary step. – whuber May 21 '13 at 16:47
  • Yes but the question, if I have understood correctly, is "Write subset-grid elevations to a file (I can't do this yet). I only need elevations (ie. Z only, not XY)" so the X,Y,Z file solution, not a grid. – gene May 21 '13 at 17:47
  • The XYZ format is a real waste here: upon removing the X and Y columns, you are left with an ASCII version of the .flt or .asc files. – whuber May 21 '13 at 18:10
  • Thanks Gene. I've just tried your suggestion (Raster > Conversion > Translate/Convert format) and I've got an XYZ ASCII file that I can manipulate for now. Indeed, I just read the file into R and calc'ed some stats on the subarea and they seem correct. Note, I realise that XYZ is very inefficient here, but for now it's a lot simpler for me to handle (baby steps to start with). – tom m May 22 '13 at 3:22
  • Then, as whuber says, it is faster to do it directly in R – gene May 22 '13 at 15:14

Here is a short R solution:

# raster data
r <- raster('filename')
# polygons, presumably from a shp file
p <- shapefile('filename.shp')
# extract values
e <- extract(r, p, df=TRUE)
# write to file
write.csv(e, 'output.txt', row.names=FALSE)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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