[Gdal-dev] ADF -> AAIGrid, "IReadBlock failed"

[Patrick Min]

Hello,

I've come up with a hacky solution to be able to
access all the data when the dataset is multitiled.

In my case, the data files are:
w001001.adf
w001000.adf
z001001.adf
z001002.adf

Below this e-mail I've included a Perl script that reads hdr.adf,
dblbnd.adf, sta.adf, and w001001x.adf, and produces some more
info than gdalinfo. Especially useful, it gives the number of tiles
in the index file, which you can use to compute where the next file
starts.

In my case w001001x.adf has info for 506520 tiles. Divide this by
HTilesPerRow (201) to get 2520, multiply by HTileYSize (4) to get
10800 rows. My HPixelSizeY is 5, so now I know the w001001.adf
has data for 10800 x 5 = 50400 pixel rows.

>From the boundary file I know that the upper left Y is at 475000,
so now I know the next file (apparently w001000.adf, strange ordering)
starts at 475000 - 50400 = 424600.

The 2nd and 3rd file have the same number of tiles as the 1st, so this
means the 3rd file starts at Y = 424600 - 50400 = 374200, and the 4th
at 374200 - 50400 = 323800. The 4th file has 170448 tiles (number
obtained by running the Perl script below as:
read_adf_info.pl z001002x.adf). All the math checks out with the
boundaries of the dataset.

To get data out of the 4th file, for example, I did the following:
- create a subdirectory "data" in the dataset directory
- move the w001* files in there
- create symbolic links:
  ln -s ./data/z001002.adf w001001.adf
  ln -s ./data/z001002x.adf w001001x.adf
  (this to make sure gdal_translate uses the 4th file)
- compute the -srcwin coordinates, taking into account the origin of
  the 4th file, which is down 3 x 50400 = 151200 from the dataset origin
- run gdal_translate

This example, more specifically: my dataset upper left is at
(13565, 475000). I wanted to extract a grid at upper left
(194990, 317925). The upper left of the 4th file is at
(13565, 323800), which, as far as gdal_translate is now concerned,
is the upper left of the whole dataset.

So the delta Y  = 323800 - 317925 = 5875, divide by cell size 5 =
5875 / 5 = row 1175.
The delta X = 194990 - 13565 = 180925 / 5 = column 36285.

The grid to extract is 600 by 600 cells, so the command becomes:

gdal_translate -oi AAIGrid -srcwin 36285 1175 600 600 hdr.adf mygrid.agr

This is a grid cell for which I had a known other height field, and
they corresponded :-)

I hope this helps people with larger datasets,
let me know if anything is not clear,

Patrick

--------------
Perl script: read_adf_info.pl
--------------
#!/usr/bin/perl

print "\nread_adf_info.pl\nwritten by Patrick Min, May 2007, v1\n\n";

$adfx_filespec = $ARGV[0];

if ($adfx_filespec eq "") {
    print "Usage: read_adf_info.pl [<adfx filespec>]\n";
    print "  if no filename is given, a default [w001001x.adf] is assumed\n\n";
    $adfx_filespec = "w001001x.adf";
}

&read_bounds_info;

&read_stats_info;

&read_header;

&read_index_info;

print "done\n\n";








sub convert_msb_int32
{
    my(@values) = @_;

    my $int_value = $values[0] * 256 * 256 * 256 + $values[1] * 256 *
256 + $values[2] * 256 + $values[3];

    return $int_value;

}  # convert_msb_int32




sub read_msb_int32
{
    my($handle) = @_;
    my $buf;

    read($handle, $buf, 4);
    my @b = unpack("CCCC", $buf);
    my $int32 = &convert_msb_int32(@b);

    return $int32;

}  # read_msb_int32



sub read_index_info
{
    #
    # construct index filename
    #
    open($index_handle, $adfx_filespec) || die "Could not open
[$adfx_filespec]\n";
    binmode $index_handle;

    #
    # read magic number
    #
    print "[$adfx_filespec] magic number [";
    read($index_handle, $magic, 8);
    @b = unpack("CCCCCCCC", $magic);
    for($i=0; $i < 8; $i++) {
	$hex = sprintf("%x", $b[$i]);
	print $hex;
	
    }
    print "]\n";

    read($index_handle, $dummy, 16);

    $filesize = 2 * &read_msb_int32($index_handle);
    print "Filesize [$filesize], ";
    $in_mb = $filesize / 1048576;
    print "  or [$in_mb MB]\n";

    read($index_handle, $dummy, 72);

    $nr_tiles = 0;
    open(LOG, "> tiles_index.log");
    while(!(eof($index_handle))) {
	
#    for($i=0; $i < $NR_TILES; $i++) {
	$tile_offset[$nr_tiles] = 2 * &read_msb_int32($index_handle);
	print LOG "Tile $nr_tiles offset [$tile_offset[$nr_tiles]]\n";

	$tile_sizes[$nr_tiles] = 2 * &read_msb_int32($index_handle);
	print LOG "Tile size [$tile_sizes[$nr_tiles]]\n";

	$nr_tiles++;
    }  # for
    close LOG;
    print "read index info for $nr_tiles tiles\n";

    close $index_handle;

}  # read_index_info




sub read_bounds_info
{
    $bounds_filespec = "dblbnd.adf";
    print "read_bounds_info($bounds_filespec)\n";

    open(BOUNDS, $bounds_filespec) || die "Could not open [$bounds_filespec]\n";
    binmode BOUNDS;

    @vars = ("LLX", "LLY", "URX", "URY");

    for($i=0; $i < 4; $i++) {

	read(BOUNDS, $buf, 8);
	@b = unpack("CCCCCCCC", $buf);

	for($j=0; $j < 8; $j++) {
	    $b2[$j] = $b[7 - $j];
	}

	$buf2 = pack("CCCCCCCC", @b2);

	$value = unpack("d", $buf2);

	print "  $vars[$i]: [$value]\n";

    }  # for

    close BOUNDS;
    print "\n";

}  # read_bounds_info



sub read_msb_double64
{
    my($handle) = @_;
    my $buf;
    my $j;
    my @b2;
    my $buf2;

    read($handle, $buf, 8);
    my @b = unpack("CCCCCCCC", $buf);

    for($j=0; $j < 8; $j++) {
	$b2[$j] = $b[7 - $j];
    }

    $buf2 = pack("CCCCCCCC", @b2);

    my $value = unpack("d", $buf2);

    return $value;

}  # read_msb_double64



sub read_stats_info
{
    $stats_filespec = "sta.adf";
    print "read_stats_info($stats_filespec)\n";

    open($stats_handle, $stats_filespec) || die "Could not open
[$stats_filespec]\n";
    binmode $stats_handle;

    @vars = ("SMin", "SMax", "SMean", "SStdDev");

    for($i=0; $i < 4; $i++) {

	$value = &read_msb_double64($stats_handle);
	print "  $vars[$i]: [$value]\n";

    }  # for

    close $stats_handle;
    print "\n";

}  # read_stats_info



sub read_header
{
    $header_filespec = "hdr.adf";
    print "read_header($header_filespec)\n";

    open($hdr_handle, $header_filespec) || die "Could not open
[$header_filespec]\n";
    binmode $hdr_handle;

    print "  HMagic [";
    for($i=0; $i < 8; $i++) {
	read($hdr_handle, $buf, 1);
	($id) = unpack("a", $buf);
	print $id;
    }
    print "]\n";

    read($hdr_handle, $dummy, 8);

    $celltype = &read_msb_int32($hdr_handle);
    print "  HCellType [$celltype]\n";

    read($hdr_handle, $dummy, 236);

    $pixel_size_x = &read_msb_double64($hdr_handle);
    print "  HPixelSizeX [$pixel_size_x]\n";
    $pixel_size_y = &read_msb_double64($hdr_handle);
    print "  HPixelSizeY [$pixel_size_y]\n";

    $unknown = &read_msb_double64($hdr_handle);
    print "  unknown1 [$unknown]\n";
    $unknown = &read_msb_double64($hdr_handle);
    print "  unknown2 [$unknown]\n";

    $tiles_per_row = &read_msb_int32($hdr_handle);
    print "  HTilesPerRow [$tiles_per_row]\n";

    $tiles_per_column = &read_msb_int32($hdr_handle);
    print "  HTilesPerColumn [$tiles_per_column]\n";

    $tile_x_size = &read_msb_int32($hdr_handle);
    print "  HTileXSize [$tile_x_size]\n";

    $unknown = &read_msb_int32($hdr_handle);
    print "  unknown [$unknown]\n";

    $tile_y_size = &read_msb_int32($hdr_handle);
    print "  HTileYSize [$tile_y_size]\n";

    close($hdr_handle);

}  # read_header