I needed an easy way to create contours to add to my topographic maps. Creating major and minor contours for each quad seemed like such a pain, as did loading individual quads into Quantum GIS, then styling them. So I wrote this script that largely automates the process, generating the minor and major (defined as 5x minor) contours from Digital Elevation Model files and combining them into one Shapefile.

I was more familiar with PHP then BASH, so I used php-cli to write this script. Writing a command line script in PHP is a bit werid, but it’s what I did. I you don’t like it, write your own script. It works well for me, but I also have PHP on my laptop for web development stuff, so your milage will vary. I store the complete set of Digital Elevation Models for NY State compressed on my laptop as BZip 2 files.

Obviously you’ll have to change paths and other minor things to make it work on your machine. Be aware you will have to install the dbase module from PECL for this script to work, as it directly reads each contour shapefile’s .dbf to convert it over to feet. I could not figure out how to do this directly with ogr2ogr, and it just seemed easier to use PHP for this task.

Output for the script is two shapefiles, the major and minor contours, as two files, consisting of all requested quads merged into major and minor shapefiles.

#!/usr/bin/php -q
<?php

array_shift($argv); // remove program name

if (!$argv || !is_numeric($argv[0]) ) echo "usage: php dbupdate.php [contour_feet] [contour1] ... [contourx]\n";

$minor = array_shift($argv);
$major = $minor*5;
$minorMeters = $minor * 0.3048;
$majorMeters = $major * 0.3048;

foreach ($argv as $quad) {

        # uncompress digital elevation model
        $dem = "/home/andy/Documents/GIS.Data/digital.elevation.model/{$quad}elu.dem";
        system("bzip2 -dkv $dem.bz2");
        
        $conDir = "/home/andy/Documents/GIS.Data/dem.contours";

        if (!file_exists($dem)) {
                echo "Digital Elevation Model for $quad NOT FOUND! \n\n";
                continue;
        }
        
        $files = glob("{$conDir}/{$quad}.{$major}ft.*");
        array_map('unlink', $files);

        $files = glob("{$conDir}/{$quad}.{$minor}ft.*");
        array_map('unlink', $files);            
                
        system("gdal_contour -a elev -i {$majorMeters} {$dem} {$conDir}/{$quad}.{$major}ft.shp");
        system("gdal_contour -a elev -i {$minorMeters} {$dem} {$conDir}/{$quad}.{$minor}ft.shp");
        
        #### convert meters to feet in minor and major contour elevations
        
        foreach (array("{$conDir}/{$quad}.{$minor}ft.dbf", "{$conDir}/{$quad}.{$major}ft.dbf") as $file) {
        
                $db = dbase_open("$file", 2);

                if ($db) {
                  $record_numbers = dbase_numrecords($db);
                  for ($i = 1; $i <= $record_numbers; $i++) {

                        // gets the old row
                        $row = dbase_get_record($db, $i);

                        // Update the date field with feet instead of meters
                        $row[1] = ceil($row[1] *3.280839);
        
                        // remove the 'deleted' entry
                        // to understand this, read the dbase man page
                        unset($row['deleted']);

                        // replace record and save
                        dbase_replace_record($db, $row, $i);  
                  }
                }

                dbase_close($db);
        }
        
        // remove uncompressed dem
        unlink($dem);           
}


// put together filename of merged file
$merged = '';
foreach ($argv as $quad) { $merged .= ".$quad"; }
$merged = 'merged'.$merged;

$quads = $argv;

$files = glob("{$conDir}/{$merged}.{$major}ft.*");
array_map('unlink', $files);

$files = glob("{$conDir}/{$merged}.{$minor}ft.*");
array_map('unlink', $files);


system("ogr2ogr {$conDir}/{$merged}.{$major}ft.shp {$conDir}/{$quads[0]}.{$major}ft.shp");
system("ogr2ogr {$conDir}/{$merged}.{$minor}ft.shp {$conDir}/{$quads[0]}.{$minor}ft.shp");

$files = glob("{$conDir}/{$quads[0]}.{$major}ft.*");
array_map('unlink', $files);

$files = glob("{$conDir}/{$quads[0]}.{$minor}ft.*");
array_map('unlink', $files);

array_shift($quads);

foreach ($quads as $quad) { 
        if (!file_exists("{$conDir}/{$quad}.{$major}ft.shp")) continue;

        system("ogr2ogr -update -append {$conDir}/{$merged}.{$major}ft.shp {$conDir}/{$quad}.{$major}ft.shp");
        
        system("ogr2ogr -update -append {$conDir}/{$merged}.{$minor}ft.shp {$conDir}/{$quad}.{$minor}ft.shp");
        
        $files = glob("{$conDir}/{$quad}.{$major}ft.*");
        array_map('unlink', $files);
        
        $files = glob("{$conDir}/{$quad}.{$minor}ft.*");
        array_map('unlink', $files);
}

?>