I need to obtain all elevation grid points from a raster DEM.
I see the steps required as:
- Read in raster DEM (Arc Info Binary Grid format).
- Create subset area using Points2One (see triangle area here: http://imgur.com/aeOPBBD).
- 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
andFortran
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: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 inR
;FOR
,DO
,WHILE
, or whatever in Fortran).R
'sraster
package supports subsetting of grids by general polygons. For rectangular subsetting, that's a basic part ofR
's array indexing capabilities.This image looks correct. The total time to read and process the dataset (before displaying it) was 1/8 second.