298 lines
7.2 KiB
Perl
Executable file
298 lines
7.2 KiB
Perl
Executable file
#!/usr/local/bin/perl
|
|
#
|
|
# dem2scene.pl -- Read in a dem data file, and output a more usable format.
|
|
#
|
|
# Written by Curtis Olson, started May 1997.
|
|
#
|
|
# Copyright (C) 1997 Curtis L. Olson - curt@infoplane.com
|
|
#
|
|
# This program is free software; you can redistribute it and/or modify
|
|
# it under the terms of the GNU General Public License as published by
|
|
# the Free Software Foundation; either version 2 of the License, or
|
|
# (at your option) any later version.
|
|
#
|
|
# This program is distributed in the hope that it will be useful,
|
|
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
# GNU General Public License for more details.
|
|
#
|
|
# You should have received a copy of the GNU General Public License
|
|
# along with this program; if not, write to the Free Software
|
|
# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
|
#
|
|
# $Id$
|
|
# (Log is kept at end of this file)
|
|
#---------------------------------------------------------------------------
|
|
|
|
|
|
use strict;
|
|
|
|
|
|
# declare variables
|
|
my($token);
|
|
my($i, $j);
|
|
my($arg);
|
|
my($res) = 1;
|
|
|
|
# "A" Record Information
|
|
my($dem_description, $dem_quadrangle);
|
|
my($dem_x1, $dem_y1, $dem_x2, $dem_y2, $dem_x3, $dem_y3, $dem_x4, $dem_y4);
|
|
my($dem_z1, $dem_z2);
|
|
my($dem_resolution, $dem_num_profiles);
|
|
|
|
# "B" Record Information
|
|
my($prof_row, $prof_col);
|
|
my($prof_num_rows, $prof_num_cols);
|
|
my($prof_x1, $prof_y1);
|
|
my($prof_data);
|
|
|
|
# set input record separator to be a space
|
|
$/ = " ";
|
|
|
|
# parse command line arguments
|
|
while ( $arg = shift(@ARGV) ) {
|
|
if ( $arg eq "-r" ) {
|
|
$res = shift(@ARGV);
|
|
if ( $res < 1 ) {
|
|
&usage();
|
|
}
|
|
} else {
|
|
&usage();
|
|
}
|
|
}
|
|
|
|
# print usage and die
|
|
sub usage {
|
|
die "Usage: $0 [ -r resval ]\n";
|
|
}
|
|
|
|
|
|
&read_a_record();
|
|
&output_scene_hdr();
|
|
|
|
$i = 0;
|
|
while ( $i < $dem_num_profiles ) {
|
|
&read_b_record();
|
|
&output_row();
|
|
|
|
$i++;
|
|
|
|
if ( $i < $dem_num_profiles ) {
|
|
# not on last record
|
|
for ( $j = 1; $j < $res; $j++ ) {
|
|
# print "skipping row\n";
|
|
&read_b_record();
|
|
}
|
|
}
|
|
}
|
|
|
|
&output_scene_close();
|
|
|
|
# read and parse DEM "A" record
|
|
sub read_a_record {
|
|
my($i);
|
|
|
|
# read initial descriptive header
|
|
while ( ($token = &next_token()) ne "_END_OF_FILE_" ) {
|
|
if ( $token !~ m/^NJ/ && $token !~ m/^NI/ ) {
|
|
$dem_description .= "$token ";
|
|
} else {
|
|
chop($dem_description);
|
|
$dem_quadrangle = $token;
|
|
last;
|
|
}
|
|
}
|
|
# print "'$dem_description' '$dem_quadrangle'\n";
|
|
|
|
# DEM level code, 3 reflects processing by DMA
|
|
&next_token();
|
|
|
|
# Pattern code, 1 indicates a regular elevation pattern
|
|
&next_token();
|
|
|
|
# Planimetric reference system code, 0 indicates geographic
|
|
# coordinate system.
|
|
&next_token();
|
|
|
|
# Zone code
|
|
&next_token();
|
|
|
|
# Map projection parameters (ignored)
|
|
for ($i = 0; $i < 15; $i++) {
|
|
&next_token();
|
|
}
|
|
|
|
# Units code, 3 represents arc-seconds as the unit of measure for
|
|
# ground planimetric coordinates throughout the file.
|
|
die "Unknown units code!\n" if ( &next_token() ne "3" );
|
|
|
|
# Units code; 2 represents meters as the unit of measure for
|
|
# elevation coordinates throughout the file.
|
|
die "Unknown units code!\n" if ( &next_token() ne "2" );
|
|
|
|
# Number (n) of sides in the polygon which defines the coverage of
|
|
# the DEM file (usually equal to 4).
|
|
die "Unknown polygon dimension!\n" if ( &next_token() ne "4" );
|
|
|
|
# Ground coordinates of bounding box in arc-seconds
|
|
$dem_x1 = &next_token();
|
|
$dem_y1 = &next_token();
|
|
|
|
$dem_x2 = &next_token();
|
|
$dem_y2 = &next_token();
|
|
|
|
$dem_x3 = &next_token();
|
|
$dem_y3 = &next_token();
|
|
|
|
$dem_x4 = &next_token();
|
|
$dem_y4 = &next_token();
|
|
|
|
# Minimum/maximum elevations in meters
|
|
$dem_z1 = &next_token();
|
|
$dem_z2 = &next_token();
|
|
|
|
# Counterclockwise angle from the primary axis of ground
|
|
# planimetric referenced to the primary axis of the DEM local
|
|
# reference system.
|
|
&next_token();
|
|
|
|
# Accuracy code; 0 indicates that a record of accuracy does not
|
|
# exist and that no record type C will follow.
|
|
# &next_token();
|
|
|
|
# DEM spacial resolution. Usually (3,3,1) (3,6,1) or (3,9,1)
|
|
# depending on latitude
|
|
$dem_resolution = &next_token();
|
|
|
|
# one dimensional arrays
|
|
&next_token();
|
|
|
|
# number of profiles
|
|
$dem_num_profiles = &next_token();
|
|
$dem_num_profiles = (($dem_num_profiles - 1) / $res) + 1;
|
|
}
|
|
|
|
|
|
# output the scene headers
|
|
sub output_scene_hdr {
|
|
my($dx, $dy, $dz);
|
|
|
|
printf("mesh %s_terrain {\n", $dem_quadrangle);
|
|
|
|
$dem_x1 =~ s/D/E/; $dem_x1 += 0.0;
|
|
$dem_y1 =~ s/D/E/; $dem_y1 += 0.0;
|
|
print " // This mesh is rooted at the following coordinates (in arc seconds)\n";
|
|
print " origin_lat = $dem_x1\n";
|
|
print " origin_lon = $dem_y1\n";
|
|
print "\n";
|
|
|
|
print " // Number of rows and columns (needed by the parser so it can create\n";
|
|
print " //the proper size structure\n";
|
|
print " rows = $dem_num_profiles\n";
|
|
print " cols = $dem_num_profiles\n"; # This isn't necessarily guaranteed
|
|
print "\n";
|
|
|
|
($dx, $dy, $dz) = $dem_resolution =~
|
|
m/(.............)(............)(............)/;
|
|
$dx *= $res;
|
|
$dy *= $res;
|
|
print " // Distance between x and y data points (in arc seconds)\n";
|
|
print " row_step = $dx\n";
|
|
print " col_step = $dy\n";
|
|
print "\n";
|
|
}
|
|
|
|
|
|
# output the scene close
|
|
sub output_scene_close {
|
|
print "\n";
|
|
print "}\n";
|
|
}
|
|
|
|
|
|
# read and parse DEM "B" record
|
|
sub read_b_record {
|
|
my($i, $j);
|
|
|
|
# row / column id of this profile
|
|
$prof_row = &next_token();
|
|
$prof_col = &next_token();
|
|
|
|
# Number of rows (elevations) and columns in this profile;
|
|
$prof_num_rows = &next_token();
|
|
$prof_num_cols = &next_token();
|
|
|
|
$prof_num_rows = (($prof_num_rows - 1) / $res) + 1;
|
|
# print "profile num rows = $prof_num_rows\n";
|
|
|
|
# Ground planimetric coordinates (arc-seconds) of the first
|
|
# elevation in the profile
|
|
$prof_x1 = &next_token();
|
|
$prof_y1 = &next_token();
|
|
|
|
# Elevation of local datum for the profile. Always zero for
|
|
# 1-degree DEM, the reference is mean sea level.
|
|
&next_token();
|
|
|
|
# Minimum and maximum elevations for the profile.
|
|
&next_token();
|
|
&next_token();
|
|
|
|
# One (usually) dimensional array ($prof_num_rows,1) of elevations
|
|
$prof_data = "";
|
|
$i = 0;
|
|
while ( $i < $prof_num_rows ) {
|
|
$prof_data .= &next_token();
|
|
$prof_data .= " ";
|
|
|
|
$i++;
|
|
|
|
if ( $i < $prof_num_rows ) {
|
|
# not on last data point
|
|
# skip the tokens to get requested resolution
|
|
for ($j = 1; $j < $res; $j++) {
|
|
# print "skipping ...\n";
|
|
&next_token();
|
|
}
|
|
}
|
|
}
|
|
chop($prof_data);
|
|
|
|
# print "$prof_data\n\n";
|
|
}
|
|
|
|
|
|
# output a row of data
|
|
sub output_row {
|
|
print " row = ($prof_data)\n";
|
|
}
|
|
|
|
|
|
# return next token from input stream
|
|
sub next_token {
|
|
my($token);
|
|
|
|
# print "in next token\n";
|
|
|
|
do {
|
|
$token = <>; chop($token);
|
|
if ( eof() ) {
|
|
$token = "_END_OF_FILE_";
|
|
}
|
|
} while ( $token eq "" );
|
|
|
|
# print "returning $token\n";
|
|
|
|
return $token;
|
|
}
|
|
|
|
while ( ($token = &next_token()) ne "_END_OF_FILE_" ) {
|
|
# print "'$token'\n";
|
|
}
|
|
|
|
|
|
#---------------------------------------------------------------------------
|
|
# $Log$
|
|
# Revision 1.1 1997/05/27 21:56:02 curt
|
|
# Initial revision (with data skipping support)
|
|
#
|