Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a high resolution digital elevation model available. The ground resolution is 1 meter and the file format is adf which is some binary ESRI grid raster format.

I have a script the outputs x/y- or lat/lon-coordinates and as a next step I have to extract the elevation of the given coordinate from the binary adf file. But that's a huge one. The file is like 1.5 GB in size and I haven't found a way yet how to access the data via a headless command line or script.

How to access adf files via script? Is there any programmes available or are there any known libraries I could use?

Or should I rather convert the adf to some ASCII text files? But that will cost a lot of disk space and performance.

I've been reading about elevation databases. How to create one from my elevation model? Is it worth a try regarding performance?

Update: To break this down, in short - How to access adf files or elevation data most efficient using a script with coordinates as parameters?

share|improve this question
1  
Are you sure its an AutoCAD format and not an Esri Grid? –  PolyGeo Jun 11 '13 at 2:44
    
Indeed, it is. I was not used to ADF and a first glance at google results yielded AutoCAD :) –  vertoe Jun 16 '13 at 8:04
    
I'm assuming that you do not have access to ArcGIS for Desktop and its Spatial Analyst extension which would be the easiest way to work with this. Have you looked at GDAL that seems to know about this format? –  PolyGeo Jun 16 '13 at 8:33
    
I have access to ArcGIS, but what I want is a script where I input x/y or lat/lon parameters and output elevation data. I tried GDAL_TRANSLATE but now I have a 8.3GB ASCII file which is as less accessable as the original ADF file. –  vertoe Jun 16 '13 at 8:39
    
I don't work with Spatial Analyst so won't be able to help personally but are you after a Python script or something else? If so, what version of ArcGIS for Desktop are you using? –  PolyGeo Jun 16 '13 at 8:55
show 1 more comment

2 Answers

One way to extract the points from adf is to open this adf file in QGIS and export it to a .xyz format and then use this file to match ur point to the point in the .xyz file

other way is to convert this adf file to a tif file using QGIS and then import this file to POSTGIS database using raster2pgsql command line tool.

this is how u import raster data (tif) to postgis Install postgresql with postgis plugin open command line and navigate to bin folder of postgresql normaly is it located here C:\Program Files\PostgreSQL\9.3\bin

STEP 1:

raster2pgsql -F -I -M -C "PATH of tif file with extension" public.your_tablename > your_output_file.sql

STEP 2 (in the same directory) excute the following

psql -U postgresql_database_Username -d Database_name -W -f your_output_file.sql

your_output_file.sql is the file generated in step 1. importing raster file to postgis is one time process n den u can use this raster table generated after importing ur tif file from the database for ur future raster queries


after this is done u can use the postgis function (ST_Value) to get the pixel value of the raster u imported

*

SELECT ST_VALUE(e.rast, ST_SetSRID(St_MakePoint(Your_Longitude, Your_Latitude ), 4326))
FROM Your_Imported_raster_table_Name e
WHERE ST_Intersects(e.rast, ST_SetSRID(St_MakePoint(Your_Longitude, Your_Latitude), 4326));

*

Hope this will help you.. i would suggest that u go for 2nd method.

share|improve this answer
add comment
up vote 0 down vote accepted

That was a quite difficult task, but now I have 2 working solutions, as described here.

First of all I had to export the binary Esri Grid file to an ASCII grid file format using the ArcToolbox > Conversion tools > From Raster > To ASCII. The result was a 2.6GB text file with a elevation data displayed in a plain text matrix with around 20k x 20k elevation values.

Now I have two different scripts I can use to access the data:

1 Single data

If I only need one single elevation data value on a certain X/Y- or LAT/LON- position I can use the following bash script:

#!/bin/bash
INPUT=$1

# Parses starting coordinates.
x11=$(awk '$1 == "xllcorner" { print $2 }' $INPUT)
y11=$(awk '$1 == "yllcorner" { print $2 }' $INPUT)

# Gets requested coordinates from args.
X=$2
Y=$3

# Calculates the line number and column.
LINE=$(bc <<< "scale=0;($Y-$y11+7)/1")
COLN=$(bc <<< "scale=0;($X-$x11+1)/1")

# Prints the Z Coordinates at X=COLN and Y=LINE
awk -v line=$LINE 'NR == line { print $0 }' $INPUT  | cut -f $COLN -d " "

Usage:

 $ ./elevation.sh dem.txt 3356387.137 5800803.818
29.58562

It takes around 7 seconds for a single elevation data to be extracted. Time consumption with this script is linear, if I call it for 100 data values, this will take around 700 seconds.

2 Multiple data

The other approach is perl script:

#!/usr/bin/perl
use strict;
use warnings;
use v5.10;
use autodie;
use List::MoreUtils qw(any);

my $data_file = shift;
my %metadata;
my @data;

open my $fh, '<', $data_file;
while (<$fh>) {
    chomp;
    my @F = split;
    if (any {$F[0] eq $_} qw(ncols nrows xllcorner yllcorner cellsize NODATA_value)) {
        $metadata{$F[0]} = $F[1];
    }
    else {
        push @data, \@F;
    }
}
close $fh;

while (@ARGV) {
    my $x = shift;
    my $y = shift;
    my $x_delta = int(($x - $metadata{xllcorner}) / $metadata{cellsize});
    my $y_delta = int(($y - $metadata{yllcorner}) / $metadata{cellsize});
    if ($x_delta < 0 or $y_delta < 0 or not defined $data[$y_delta][$x_delta]) {
        say $metadata{NODATA_value};
    }
    else {
        say $data[$y_delta][$x_delta];
    }
}

Usage with multiple coordinates:

 $ ./elevation.pl dem.txt  0 0 3356385.137 5800799.818 3356387.137 5800803.818
-9999
31.11266
29.58562

It takes around 2 minutes to compute which is much longer than the bash script. But the time consumption of this script is not connected to the number of arguments. To access 100 elevation values it still takes only 2 minutes.

Downside of this script it needs alot of main memory, around 10x of the DEM file size, in my case: 26GB.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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