google-elev-api.pl to HTML.

index -|- end

Generated: Mon Aug 29 19:34:40 2016 from google-elev-api.pl 2016/07/15 34.6 KB. text copy

#!/usr/bin/perl -w
# NAME: google-elev-api.pl
# AIM: SPECIFIC - Read a downloaded Google API elevation file
use strict;
use warnings;
use File::Basename;  # split path ($name,$dir,$ext) = fileparse($file [, qr/\.[^.]*/] )
use Math::Trig;
use Cwd;
my $os = $^O;
my $perl_dir = '/home/geoff/bin';
my $PATH_SEP = '/';
my $temp_dir = '/tmp';
if ($os =~ /win/i) {
    $perl_dir = 'C:\GTools\perl';
    $temp_dir = $perl_dir;
    $PATH_SEP = "\\";
}
unshift(@INC, $perl_dir);
require 'lib_utils.pl' or die "Unable to load 'lib_utils.pl' Check paths in \@INC...\n";
# log file stuff
our ($LF);
my $pgmname = $0;
if ($pgmname =~ /(\\|\/)/) {
    my @tmpsp = split(/(\\|\/)/,$pgmname);
    $pgmname = $tmpsp[-1];
}
my $outfile = $temp_dir.$PATH_SEP."temp.$pgmname.txt";
open_log($outfile);

# user variables
my $VERS = "0.0.5 2015-01-09";
my $load_log = 0;
my $in_file = '';
my $verbosity = 0;
my $out_file = '';
my $add_blue_cont = 0;
my $add_green_cont = 0;
my $use_bl_of_tile = 0; # def use x,y as TL, x+1,y+1 BR
my $add_contours = 0;

####################################
###############
my $def_lat = 27.687822920;    # 27.687822920
my $def_lon = 86.728243490;    # 86.72824349;
#my $def_zoom = 14;
##my $def_lat = 27.6730;
##my $def_lon = 86.7171;
my $def_zoom = 12;
####################################

my $cont_cnt = 5;   # 20;  # 10;
my $url_base = "http://a.tile.openstreetmap.org";
my $xg_out = $temp_dir.$PATH_SEP."tempslip2.xg";

# Theo's values
# my $def_samples = 30;
# my $def_tile_X = 12139;
# my $def_tile_Y = 6879;
my $UL_lat = 27.7613298745; # 27.6835280837878; # 27.702983735525837;
my $UL_lon = 86.6601562500; # 86.66015625;      # 86.72607421875;
my $LR_lat = 27.6835280838; # 27.6056708264655; # 27.683528083787767;
my $LR_lon = 86.7480468750; # 86.748046875;     # 86.748046875;

#static const char * xgraph_colors[] = 
#    {"black", "white", "red", "blue", "green", "violet", // 0,  ..., 5
#     "orange", "yellow", "pink", "cyan", "#A2B5CD",      // 6,  ..., 10
#     "#6C7B8B", "#FF00FF", "#00CDCD", "navy", "gold"     // 11, ..., 15
#    };

#my @xg_colors = qw( black white red blue green violet 
#    orange yellow pink cyan #A2B5CD 
#    #6C7B8B" #FF00FF #00CDCD navy gold );

# ### DEBUG ###
my $debug_on = 1;
my $def_file = 'C:\Users\user\Downloads\elevations_28_87-z14-t1-30x30.txt';
# Slippy Tile: http://b.tile.openstreetmap.org/12/3034/1719.png
my $def_dir = 'F:\FGx\fgx.github.io\sandbox\flightpaths\vnlk';
my $def_file2 = $def_dir.'\elevations_27.6730_86.7171_z12_t5_250x250_vnlk-both-paths.txt';
my $def_file3 = $def_dir.'\elevations_27.7111_86.7123_z12_t4_500x500_vnlk.txt';
my $def_file4 = $def_dir.'\elevations_27.7224_86.7041_z12_t3_100x100_vnlk.txt';
my $def_file5 = 'F:\FGx\fgx.github.io\sandbox\flightpaths\vhsk\elevations_22.4366_114.0804_z14_t3_200x200_.txt';

my $def_tiles = 5;

my $flightPathFile = $def_dir.'\VNLK-01-cooked.csv';    # red
my $flightPathFile2 = $def_dir.'\VNLK-02-cooked.csv';   # green
my $flightPathFile3 = $def_dir.'\6-25-2016-1-cooked.csv'; # purple

### program variables
my @warnings = ();
my $cwd = cwd();
my $xg = "# gen $pgmname, on ".lu_get_YYYYMMDD_hhmmss_UTC(time())." UTC\n";
my $g_rmm_tile;
my ($g_imin,$g_imax);

my $M2FT = 3.28084;

sub VERB1() { return $verbosity >= 1; }
sub VERB2() { return $verbosity >= 2; }
sub VERB5() { return $verbosity >= 5; }
sub VERB9() { return $verbosity >= 9; }

sub show_warnings($) {
    my ($val) = @_;
    if (@warnings) {
        prt( "\nGot ".scalar @warnings." WARNINGS...\n" );
        foreach my $itm (@warnings) {
           prt("$itm\n");
        }
        prt("\n");
    } else {
        prt( "\nNo warnings issued.\n\n" ) if (VERB9());
    }
}

sub pgm_exit($$) {
    my ($val,$msg) = @_;
    if (length($msg)) {
        $msg .= "\n" if (!($msg =~ /\n$/));
        prt($msg);
    }
    show_warnings($val);
    close_log($outfile,$load_log);
    exit($val);
}


sub prtw($) {
   my ($tx) = shift;
   $tx =~ s/\n$//;
   prt("$tx\n");
   push(@warnings,$tx);
}

###########################################################

sub getTileNumber($$$) {
    my ($lat,$lon,$zoom) = @_;
    my $z2 = 2 ** $zoom;
    my $xtile = int( ($lon+180)/360 * $z2 ) ;
    my $ytile = int( (1 - log(tan(deg2rad($lat)) + sec(deg2rad($lat)))/pi)/2 * $z2 ) ;
    return ($xtile, $ytile);
}

sub getLonLat {
   my ($xtile, $ytile, $zoom) = @_;
   my $n = 2 ** $zoom;
   my $lon_deg = $xtile / $n * 360.0 - 180.0;
   my $lat_deg = rad2deg(atan(sinh(pi * (1 - 2 * $ytile / $n))));
   return ($lon_deg, $lat_deg);
}


##########################################################
sub get_double_stg($) {
    my $d = shift;
    my $cd = sprintf("%4.10f",$d);
    return $cd;
}

sub my_round($) {
    my $d = shift;
    my $v = int( $d + 0.5 );
    $v = int( $d - 0.5 ) if ($d < 0);
    return $v;
}

##########################################################
sub get_minmax_ref() {
    my @a = (200,200,-200,-200);
    return \@a;
}
sub set_minmax_ref($$$) {
    my ($rmm,$lat,$lon) = @_;
    ${$rmm}[0] = $lon if ($lon < ${$rmm}[0]);
    ${$rmm}[1] = $lat if ($lat < ${$rmm}[1]);
    ${$rmm}[2] = $lon if ($lon > ${$rmm}[2]);
    ${$rmm}[3] = $lat if ($lat > ${$rmm}[3]);
}

sub expand_minmax_ref($$) {
    my ($rmm,$factor) = @_;
    return if ($factor <= 0);
    my $min_lon = ${$rmm}[0];
    my $min_lat = ${$rmm}[1];
    my $max_lon = ${$rmm}[2];
    my $max_lat = ${$rmm}[3];
    my $lon_dif = $max_lon - $min_lon;
    my $lat_dif = $max_lat - $min_lat;
    # do we have positive difference
    return if ( !(($lon_dif > 0)&&($lat_dif > 0)) );
    # calculate the smudge factors - in degrees
    my $lon_sm = (($factor / 100) * $lon_dif) / 2;
    my $lat_sm = (($factor / 100) * $lat_dif) / 2;
    # apply the smudge factors
    $min_lon -= $lon_sm;
    $min_lat -= $lat_sm;
    $max_lon += $lon_sm;
    $max_lat += $lat_sm;
    # update the reference
    set_minmax_ref($rmm,$min_lat,$min_lon);
    set_minmax_ref($rmm,$max_lat,$max_lon);
}

sub get_minmax_ref_xg($$) {
    my ($rmm,$color) = @_;
    my $min_lon = ${$rmm}[0];
    my $min_lat = ${$rmm}[1];
    my $max_lon = ${$rmm}[2];
    my $max_lat = ${$rmm}[3];
    my $xg = "# BBOX: $min_lon,$min_lat,$max_lon,$max_lat\n";
    $xg .= "color $color\n";
    $xg .= "$min_lon $min_lat\n";
    $xg .= "$min_lon $max_lat\n";
    $xg .= "$max_lon $max_lat\n";
    $xg .= "$max_lon $min_lat\n";
    $xg .= "$min_lon $min_lat\n";
    $xg .= "NEXT\n";
    return $xg;
}

sub get_latlon_xy($$$$) {
    my ($rmm,$x,$y,$sqr) = @_;
    my $tl_lon = ${$rmm}[0]; # $min_lon
    my $tl_lat = ${$rmm}[3]; # $max_lat
    my $br_lon = ${$rmm}[2]; # $max_lon
    my $br_lat = ${$rmm}[1]; # $min_lat
    my $lon = (($x / $sqr) * ($br_lon - $tl_lon)) + $tl_lon; # add $min_lon
    my $lat = (($y / $sqr) * ($tl_lat - $br_lat)) + $br_lat; # add $min_lat
    return ($lat,$lon);
}

###########################################################################
sub get_theo_xg() {
    my $xg = "# Theo lon lat UL $UL_lon $UL_lat LR $LR_lon $LR_lat\n";
    my $t_ctr_lon = ($UL_lon + $LR_lon) / 2;
    my $t_ctr_lat = ($UL_lat + $LR_lat) / 2;
    $xg .= "color red\n";
    $xg .= "$UL_lon $UL_lat\n";
    $xg .= "$LR_lon $LR_lat\n";
    $xg .= "NEXT\n";
    $xg .= "color green\n";
    $xg .= "$t_ctr_lon $t_ctr_lat\n";
    $xg .= "NEXT\n";
    $xg .= "anno $t_ctr_lon $t_ctr_lat Tile Ctr\n";
    return $xg;
}

sub get_slippy_tiles($$$$) {
    my ($lat,$lon,$z,$t) = @_;
    my ($x,$y) = getTileNumber($lat,$lon,$z);
    my $dir = "/$z/$x/$y.png";
    my $tc = int($t / 2);
    my ($xi,$yi);
    # get TOP LEFT
    my $x1 = $x - $tc;
    my $y1 = $y - $tc;
    # get BOT RIGHT
    my $x2 = $x + $tc;
    my $y2 = $y + $tc;
    prt("PNG: $url_base$dir - $tc\n");
    $xg .= "# PNG: $url_base$dir\n";
    
    my ($tl_lon, $tl_lat);
    my ($bl_lon, $bl_lat);
    my ($tr_lon, $tr_lat);
    my ($br_lon, $br_lat);
    for (; $x1 <= $x2; $x1++) {
        $y1 = $y - $tc;
        for (; $y1 <= $y2; $y1++) {
            ($tl_lon, $tl_lat) = getLonLat($x1+0,$y1+0,$z);
            ($bl_lon, $bl_lat) = getLonLat($x1+0,$y1+1,$z);
            ($tr_lon, $tr_lat) = getLonLat($x1+1,$y1+0,$z);
            ($br_lon, $br_lat) = getLonLat($x1+1,$y1+1,$z);
            my $rmm = get_minmax_ref();
            set_minmax_ref($rmm,$tl_lat,$tl_lon);
            set_minmax_ref($rmm,$bl_lat,$bl_lon);
            set_minmax_ref($rmm,$tr_lat,$tr_lon);
            set_minmax_ref($rmm,$br_lat,$br_lon);
            $xg .= get_minmax_ref_xg($rmm,'blue');
        }
    }

    $x1 = $x - $tc;
    $y1 = $y - $tc;
    $x2 = $x + $tc;
    $y2 = $y + $tc;
    ($tl_lon, $tl_lat) = getLonLat($x-$tc,$y-$tc,$z);
    ($bl_lon, $bl_lat) = getLonLat($x-$tc,$y+$tc+1,$z);
    ($tr_lon, $tr_lat) = getLonLat($x+$tc+1,$y+$tc,$z);
    ($br_lon, $br_lat) = getLonLat($x+$tc+1,$y-$tc+1,$z);
    my $ctr_lat = ($tl_lat + $br_lat) / 2;
    my $ctr_lon = ($tl_lon + $br_lon) / 2;
    my ($cx, $cy) = getTileNumber($ctr_lat,$ctr_lon,$z);
    $xg .= "anno $ctr_lon $ctr_lat Cntr Elevs $cx,$cy\n";
    $xg .= "anno $lon $lat Tile Ctr $x,$y\n";
    my $rmmt = get_minmax_ref();
    set_minmax_ref($rmmt,$tl_lat,$tl_lon);
    set_minmax_ref($rmmt,$bl_lat,$bl_lon);
    set_minmax_ref($rmmt,$tr_lat,$tr_lon);
    set_minmax_ref($rmmt,$br_lat,$br_lon);
    $g_rmm_tile = $rmmt;
    $xg .= get_minmax_ref_xg($rmmt,'gray');

}

sub get_slippy_tile($$$) {
    my ($lat,$lon,$z) = @_;

    my ($x,$y) = getTileNumber($lat,$lon,$z);
    my $dir = "/$z/$x/$y.png";
    prt("PNG: $url_base$dir\n");
    $xg .= "# PNG: $url_base$dir\n";
    
    my ($tl_lon, $tl_lat);
    my ($bl_lon, $bl_lat);
    my ($tr_lon, $tr_lat);
    my ($br_lon, $br_lat);

    if ($use_bl_of_tile) {
        ($tl_lon, $tl_lat) = getLonLat($x+0,$y+1,$z);
        ($bl_lon, $bl_lat) = getLonLat($x+0,$y+0,$z);
        ($tr_lon, $tr_lat) = getLonLat($x+1,$y+1,$z);
        ($br_lon, $br_lat) = getLonLat($x+1,$y+0,$z);
    } else {
        ($tl_lon, $tl_lat) = getLonLat($x+0,$y+0,$z);
        ($bl_lon, $bl_lat) = getLonLat($x+0,$y+1,$z);
        ($tr_lon, $tr_lat) = getLonLat($x+1,$y+0,$z);
        ($br_lon, $br_lat) = getLonLat($x+1,$y+1,$z);
    }

    my $ctr_lat = ($tl_lat + $br_lat) / 2;
    my $ctr_lon = ($tl_lon + $br_lon) / 2;
    my ($cx, $cy) = getTileNumber($ctr_lat,$ctr_lon,$z);

    my $rmmt = get_minmax_ref();
    set_minmax_ref($rmmt,$tl_lat,$tl_lon);
    set_minmax_ref($rmmt,$bl_lat,$bl_lon);
    set_minmax_ref($rmmt,$tr_lat,$tr_lon);
    set_minmax_ref($rmmt,$br_lat,$br_lon);
    $g_rmm_tile = $rmmt;
    $xg .= get_minmax_ref_xg($rmmt,'gray');

    $xg .= get_theo_xg();

}

sub get_track($$) {
    my ($inf,$color) = @_;
    if (! open INF, "<$inf") {
        pgm_exit(1,"ERROR: Unable to open file [$inf]\n"); 
    }
    my @lines = <INF>;
    close INF;
    my $lncnt = scalar @lines;
    prt("Processing $lncnt lines, from [$inf]...\n");
    my ($line,@arr,$lon,$lat,$alt,$lon1,$lat1,$alt1,$cnt,$half);
    $cnt = 0;
    my $xg = "color $color\n";
    $half = int($lncnt / 2);
    foreach $line (@lines) {
        $cnt++;
        chomp $line;
        @arr = split(/,/,$line);
        next if ($arr[0] eq 'lon');
        $lon = $arr[0];
        $lat = $arr[1];
        $alt = $arr[2];
        $xg .= "$lon $lat ; $alt\n";
        if ($cnt == $half) {
            $xg .= "anno $lon $lat Track $color\n";
        }
    }
    $xg .= "NEXT\n";
    return $xg;
}


sub write_xg() {
    $xg .= get_VNLK_xg();

    $xg = get_track($flightPathFile,'red').$xg;
    $xg = get_track($flightPathFile2,'green').$xg;
    $xg = get_track($flightPathFile3,'violet').$xg;

    rename_2_old_bak($xg_out);
    write2file($xg,$xg_out);
    prt("Written xg to $xg_out\n");

}

sub get_contours($$$$$);


# expect elev _ lat     _ lon     _ zoom _ tiles _ sqrxsqr _ other...
# 0             1         2         3      4       5         6
# elevations  _ 27.6730 _ 86.7171 _ z12  _ t5    _ 250x250 _ vnlk-both-paths.txt
sub parse_filename($$$$$$) {
    my ($inf,$rlat,$rlon,$rzoom,$rtiles,$rsqr) = @_;
    my ($n,$d,$e) = fileparse($inf, qr/\.[^.]*/ );
    my @arr = split("_",$n);
    my $cnt = scalar @arr;
    if ($cnt < 6) {
        pgm_exit(1,"File name '$n' split to $cnt! Expect min 6...\n");
    }
    my $lat = $arr[1];
    my $lon = $arr[2];
    my $zm  = $arr[3];
    my $z = 0;
    if ($zm =~ /^z(\d+)$/) {
        $z = $1;
    } else {
        pgm_exit(1,"File name '$n' split $cnt, but 3 '$zm' not zoom!\n");
    }
    my $tcnt = $arr[4];
    my $t = 0;
    if ($tcnt =~ /^t(\d+)$/) {
        $t = $1;
    } else {
        pgm_exit(1,"File name '$n' split $cnt, but 4 '$tcnt' not tile count!\n");
    }
    my $sqrxsqr = $arr[5];
    my @a = split("x",$sqrxsqr);
    $tcnt = scalar @a;
    my $sqr = $a[0];
    if ($tcnt != 2) {
        pgm_exit(1,"File name '$n' split $cnt, but 5 '$sqrxsqr' not sqr x sqr!\n");
    }
    if ($sqr != $a[1]) {
        pgm_exit(1,"File name '$n' split $cnt, but 5 '$sqr' not = $a[1]!\n");

    }
    # this is the 'center' of a tile in the tile array
    get_slippy_tiles($lat,$lon,$z,$t);
    ${$rlat} = $lat;
    ${$rlon} = $lon;
    ${$rzoom} = $z;
    ${$rtiles} = $t;
    ${$rsqr} = $sqr;

    prt("Name '$n' = $lat,$lon,$z,$t,$sqr which seems ok...\n");
    # write_xg();
    # pgm_exit(1,"TEMP EXIT\n");
}

sub process_in_file($) {
    my ($inf) = @_;
    if (! open INF, "<$inf") {
        pgm_exit(1,"ERROR: Unable to open file [$inf]\n"); 
    }
    my @lines = <INF>;
    close INF;
    my $lncnt = scalar @lines;
    my ($flat,$flon,$fzoom,$ftiles,$fsqr);
    parse_filename($inf,\$flat,\$flon,\$fzoom,\$ftiles,\$fsqr);
    prt("Processing $lncnt lines, from [$inf]...\n");
    my ($line,$cnt,$lnn,@arr,$sqr,$alt,$x,$y,$pos,$rng,$step,$i,$j,$tx,$ty);
    my ($lat,$lon,$feet,$nxt,$j2,$tmp,$ra);
    my $min_alt =  999999;
    my $max_alt = -999999;
    $lnn = 0;
    my @steps = ();
    my @arrxy = ();
    #foreach $line (@lines) {
    $line = $lines[0];
    chomp $line;
    $lnn++;
    @arr = split(",",$line);
    $cnt = scalar @arr;
    $sqr = sqrt($cnt);
    prt("$cnt elevations - $sqr x $sqr (=$fsqr?)\n");
    $x = 0;
    $y = 0;
    my $mean_alt = 0;
    # First run through the array of ELEVATIONS
    # and setup an array of arrays for x,y addressing
    # Also get min and max elevations
    for ($j = 0; $j < $cnt; $j++) {
        $alt = $arr[$j];
        $mean_alt += int($alt);
        if ($alt < $min_alt) {
            $min_alt = $alt;
            $g_imin = $j;
            $ty = int($j / $sqr);
            $tx = $j - ($ty * $sqr);
            if (($x == $tx)&&($y == $ty)) {
            } else {
                pgm_exit(1,"Calc FAILED [$x != $tx] or [$y != $ty] on $j\n");
            }
        }
        if ($alt > $max_alt) {
            $max_alt = $alt;
            $g_imax = $j;
            $ty = int($j / $sqr);
            $tx = $j - ($ty * $sqr);
            if (($x == $tx)&&($y == $ty)) {
            } else {
                pgm_exit(1,"Calc FAILED [$x != $tx] or [$y != $ty] on $j\n");
            }
        }

        $arrxy[$x][$y] = $alt;

        # bump the offsets
        ###################################
        $x++;
        if ($x >= $sqr) {
            $x = 0;
            $y++;
        }
        ###################################
    }
    $mean_alt = int($mean_alt / $cnt);
    $rng = $max_alt - $min_alt;
    $step = $rng / $cont_cnt;
    for ($i = 0; $i <= $cont_cnt; $i++) {
        $alt = $min_alt + ($i * $step);
        $steps[$i] = my_round($alt);
    }
    prt("Min $min_alt, Max $max_alt, Mean $mean_alt, Diff $rng, step $step\n");
    my $min_feet = my_round($min_alt * $M2FT); # get the feet
    my $max_feet = my_round($max_alt * $M2FT);
    prt("Steps: $min_alt m $min_feet ft: ");
    for ($i = 0; $i <= $cont_cnt; $i++) {
        $alt = $steps[$i];
        $feet = my_round( $alt * $M2FT );
        prt("$feet ");
    }
    prt(": $max_alt m. $max_feet ft.");
    prt("\n");

    $xg .= "# Add min/max\n";
    $ty = int($g_imin / $sqr);
    $tx = $g_imin - ($ty * $sqr);
    $x = $tx;
    $y = $ty;
    $pos = ($y * $sqr) + $x;
    $alt = $arr[$pos];
    ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x,$y,$sqr);
    prt("At $x,$y($pos) = $alt, $lat, $lon\n");
    $xg .= "anno $lon $lat Min $alt ($x,$y)\n";
    $xg .= "color yellow\n";
    $xg .= "$lon $lat ; $alt ($x $y)\n";
    $xg .= "NEXT\n";

    $ty = int($g_imax / $sqr);
    $tx = $g_imax - ($ty * $sqr);
    $x = $tx;
    $y = $ty;
    $pos = ($y * $sqr) + $x;
    $alt = $arr[$pos];
    ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x,$y,$sqr);
    prt("At $x,$y($pos) = $alt, $lat, $lon\n");
    $xg .= "anno $lon $lat Max $alt ($x,$y)\n";
    $xg .= "color yellow\n";
    $xg .= "$lon $lat ; $alt ($x $y)\n";
    $xg .= "NEXT\n";

    # try to add contour for MEAN altitude
    my @points = ();
    for ($j = 0; $j < $cnt; $j++) {
        $alt = $arr[$j];
        $j2 = $j + 1;
        if ($j2 < $cnt) {
            $nxt = $arr[$j2];
            if ($alt > $nxt) {
                $tmp = $alt;
                $alt = $nxt;
                $nxt = $tmp;
            }
            if (($mean_alt >= $alt) && ($mean_alt <= $nxt)) {
                push(@points,[$j,$j2]);
            }
        }
    }
    $nxt = scalar @points;
    prt("Got $nxt points...\n");
    if ($nxt && $add_green_cont) {
        $xg .= "# got $nxt points...\n";
        $xg .= "color green\n";
        foreach $ra (@points) {
            $j = ${$ra}[0];
            $j2 = ${$ra}[1];
            $alt = $arr[$j];
            $y = int($j / $sqr);
            $x = $j - ($y * $sqr);
            $tmp = $arrxy[$x][$y];
            ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x,$y,$sqr);
            # prt("At $x,$y($pos) = $alt, $lat, $lon\n");
            $xg .= "$lon $lat ; $alt ($x $y) $tmp\n";
            $ty = int($j2 / $sqr);
            $tx = $j2 - ($y * $sqr);
            $tmp = $arrxy[$tx][$ty];
            ($lat,$lon) = get_latlon_xy($g_rmm_tile,$tx,$ty,$sqr);
            if (($x == $tx)||($y == $ty)) {
                $xg .= "$lon $lat ; $tmp ($tx $ty)\n";
            }
            $xg .= "NEXT\n";
        }
    }

    # try stepping vertically, using x,y offsets
    my @pts2 = ();
    for ($x = 0; $x < $sqr; $x++) { # for each col
        for ($y = 0; $y < $sqr; $y++) { # go down the rows
            $ty = $y + 1;
            $alt = $arrxy[$x][$y];
            if ($ty < $sqr) {
                $nxt = $arrxy[$x][$ty];
                if ($alt > $nxt) {
                    $tmp = $alt;
                    $alt = $nxt;
                    $nxt = $tmp;
                }
                if (($mean_alt >= $alt) && ($mean_alt <= $nxt)) {
                    $j = $y * $sqr + $x;
                    $j2 = ($ty * $sqr) + $x;
                    push(@pts2,[$j,$j2]);
                }
            }
        }
    }
    $nxt = scalar @pts2;
    prt("Got $nxt points...\n");
    if ($nxt && $add_blue_cont) {
        $xg .= "# got $nxt points...\n";
        $xg .= "color blue\n";
        foreach $ra (@pts2) {
            $j = ${$ra}[0];
            $j2 = ${$ra}[1];
            $alt = $arr[$j];
            $y = int($j / $sqr);
            $x = $j - ($y * $sqr);
            $tmp = $arrxy[$x][$y];
            ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x,$y,$sqr);
            # prt("At $x,$y($pos) = $alt, $lat, $lon\n");
            $xg .= "$lon $lat ; $alt ($x $y) $tmp\n";
            $ty = int($j2 / $sqr);
            $tx = $j2 - ($ty * $sqr);
            $tmp = $arrxy[$tx][$ty];
            ($lat,$lon) = get_latlon_xy($g_rmm_tile,$tx,$ty,$sqr);
            if (($x == $tx)||($y == $ty)) {
                $xg .= "$lon $lat ; $tmp ($tx $ty)\n";
            }
            $xg .= "NEXT\n";
        }
    }

    if ($add_contours) {
        my $color = 3;
        my $ptcnt = 0;
        for ($i = 0; $i <= $cont_cnt; $i++) {
            $alt = $steps[$i];
            $ptcnt += get_contours( \@arrxy, $sqr, $alt, $color, 1 );
            $color++;
        }
        prt("Added $ptcnt contour points...\n");
        $alt = $mean_alt;
        # get_contours( \@arrxy, $sqr, $alt, "white", 1 );
    }

}

sub get_VNLK_xg() {
    my $txt = <<EOF;
# icao="VNLK", name:"Lukla", lon:086.72824349, lat:27.68782292, alt:9100, rwys:1, ids:"24/06"
anno 086.72824349 27.68782292 VNLK Twr: 
# Gen VNLK CIRCUIT AI 3 gs, alt 1000 ft.(305m), Km dist 5.816, 141 secs at 41 mps, 80 kts ias, -424 fpm
# NO bbox!
# begin runway description
anno 086.72545143 27.68513766 06
anno 086.72914112 27.68807225 24
color blue
anno 86.7272962502131 27.6866049673658 len:1601 ft
086.72545143 27.68513766
086.72914112 27.68807225
NEXT
color red
86.7255523043977 27.6850371581307
86.7253505554176 27.6852381617944
86.7290402427208 27.6881727517518
86.7292419970945 27.6879717481733
86.7255523043977 27.6850371581307
NEXT
# end runway description

EOF
    return $txt;
}

###################################
### How to arrange these points into line

sub ras_have_same_side($$) {
    my ($ra1,$ra2) = @_;
    my $x1 = ${$ra1}[0];
    my $y1 = ${$ra1}[1];
    my $x2 = ${$ra2}[0];
    my $y2 = ${$ra2}[1];
    if ($x1 == $x2) {
        return 1 if ($y1 == $y2);
        return 1 if (abs($y1 - $y2) == 1);
    } elsif ($y1 == $y2) {
        return 1 if (abs($x1 - $x2) == 1);
    }
    #prt("Test ($x1,$y1) - ($x2,$y2)\n");
    return 0;
}



sub sort_matches($) {
    my $rma = shift; # =  \@matches
    my ($ra1,$ra2,$x1,$x2,$y1,$y2,$i,$j);
    #my ($x11,$x21,$y11,$y21);
    my $max = scalar @{$rma};
    return if ($max < 3);
    prt("Sort $max matches into lines...\n");
    $ra1 = ${$rma}[0];
    $x1 = ${$ra1}[0];
    $y1 = ${$ra1}[1];
    my @nmatches = ();
    my %added = ();
    my $cnt = 0;
    for ($i = 0; $i < $max; $i++) {
        next if (defined $added{$i});
        $ra1 = ${$rma}[$i];
        # point 1
        $x1 = ${$ra1}[0];
        $y1 = ${$ra1}[1];
        # expand to square
        #$x11 = $x1 + 1;
        #$y11 = $y1 + 1;
        # seek a match == share a side
        for ($j = 0; $j < $max; $j++) {
            next if ($i == $j);
            next if (defined $added{$j});
            $ra2 = ${$rma}[$j];
            $x2 = ${$ra2}[0];
            $y2 = ${$ra2}[1];
            #$x21 = $x2 + 1;
            #$y21 = $y2 + 1;
            if (ras_have_same_side($ra1,$ra2)) {
                prt("ras_have_same_side ($x1,$y1) - ($x2,$y2) - $i and $j\n") if (VERB5());
                if (! defined $added{$i} ) {
                    push(@nmatches,$ra1);
                    $cnt++;
                }
                push(@nmatches,$ra2);
                $cnt++;
                $added{$i} = 1;
                $added{$j} = 1;
            }
        }
    }
    for ($i = 0; $i < $max; $i++) {
        next if (defined $added{$i});
        $ra1 = ${$rma}[$i];
        push(@nmatches,$ra1);
    }
    prt("Got $cnt of $max matched into lines...\n");
    # $load_log = 1;
    @{$rma} = @nmatches;
}

sub get_contours($$$$$) {
    my ($raxy,$sqr,$alt,$color,$flag) = @_;
    my $max = $sqr - 1;
    return if ($max < 5);
    my ($x1,$x2,$y1,$y2,$tx,$ty);
    my ($alt1,$alt2,$alt3,$alt4,$tmp);
    my ($lat,$lon,$match,$xm,$ym);
    my ($min_alt,$max_alt,$av_alt,$ra,$adif,$min_dif);
    my @matches = ();
    $xg .= "color $color\n";
    $match = 0;
    for ($x1 = 0; $x1 < $max; $x1++) { # for across the cols
        $x2 = $x1 + 1;
        for ($y1 = 0; $y1 < $max; $y1++) { # go down the rows
            $y2 = $y1 + 1;
            $alt1 = ${$raxy}[$x1][$y1]; # LT
            $alt2 = ${$raxy}[$x2][$y1]; # RT
            $alt3 = ${$raxy}[$x2][$y2]; # RB
            $alt4 = ${$raxy}[$x1][$y2]; # LB
            ### get min/max of grid
            $min_alt = $alt1;
            $max_alt = $alt1;
            $min_alt = $alt2 if ($alt2 < $min_alt);
            $max_alt = $alt2 if ($alt2 > $max_alt);
            $min_alt = $alt3 if ($alt3 < $min_alt);
            $max_alt = $alt3 if ($alt3 > $max_alt);
            $min_alt = $alt4 if ($alt4 < $min_alt);
            $max_alt = $alt4 if ($alt4 > $max_alt);
            ### does altitude fit in this grid
            if (($alt >= $min_alt) && ($alt <= $max_alt)) {
                push(@matches,[$x1,$y1]);
                $match++;
            }
        }
    }
    ###################################
    ### How to arrange these points into line
    sort_matches(\@matches);
    $max = scalar @matches;
    prt("# Found $max pts, with mean $alt\n");
    $xg .= "# Found $max pts, with mean $alt\n";
    my $f = 0.1;
    my ($i,$i2,$ra2);
    my $open = 0;
    for ($i = 0; $i < $max; $i++) {
        $ra = $matches[$i];
        $i2 = $i + 1;
        $x1 = ${$ra}[0];
        $y1 = ${$ra}[1];
        $x2 = $x1 + 1;
        $y2 = $y1 + 1;
        $alt1 = ${$raxy}[$x1][$y1]; # LT
        $alt2 = ${$raxy}[$x2][$y1]; # RT
        $alt3 = ${$raxy}[$x2][$y2]; # RB
        $alt4 = ${$raxy}[$x1][$y2]; # LB
        ###########################
        ### Diffs
        $adif = abs($alt - $alt1);
        $xm = $x1 + $f;
        $ym = $y1 + $f;
        $min_dif = $adif;

        $adif = abs($alt - $alt2);
        if ($adif < $min_dif) {
            $min_dif = $adif;
            $xm = $x2 - $f;
            $ym = $y1 + $f;
        }
        $adif = abs($alt - $alt3);
        if ($adif < $min_dif) {
            $min_dif = $adif;
            $xm = $x2 - $f;
            $ym = $y2 - $f;
        }
        $adif = abs($alt - $alt4);
        if ($adif < $min_dif) {
            $min_dif = $adif;
            $xm = $x1 + $f;
            $ym = $y2 - $f;
        }
        ################################
        ### Paint the point
        ($lat,$lon) = get_latlon_xy($g_rmm_tile,$xm,$ym,$sqr);
        $xg .= "$lon $lat ; $alt1 $alt2 $alt3 $alt4 ($xm,$ym)\n";
        if ($i2 < $max) {
            $ra2 = $matches[$i2];
            if (ras_have_same_side($ra,$ra2)) {
                $open = 1;
            } else {
                $xg .= "NEXT\n";
                $open = 0;
            }
        } else {
            $xg .= "NEXT\n";
            $open = 0;
        }
        if ($flag & 0x80) {
            ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x1,$y1,$sqr);
            $xg .= "$lon $lat ; $alt1\n";
            ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x2,$y1,$sqr);
            $xg .= "$lon $lat ; $alt2\n";
            ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x2,$y2,$sqr);
            $xg .= "$lon $lat ; $alt3\n";
            ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x1,$y2,$sqr);
            $xg .= "$lon $lat ; $alt4\n";
            ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x1,$y1,$sqr);
            $xg .= "$lon $lat ; $alt1\n";
            $xg .= "NEXT\n";
        }
    }
    $xg .= "NEXT\n" if ($open);
    return $max;
}

sub get_contours_old($$$$$) {
    my ($raxy,$sqr,$alt,$color,$flag) = @_;
    my ($x1,$x2,$y1,$y2,$tx,$ty);
    my ($alt1,$alt2,$alt3,$alt4,$tmp);
    my ($altl,$alth,$fnd);
    my ($lat,$lon,$match,$xm,$ym);
    my $max = $sqr - 1;
    my @matches = ();
    $xg .= "color $color\n";
    $fnd = 0;
    for ($x1 = 0; $x1 < $max; $x1++) { # for each col
        $x2 = $x1 + 1;
        for ($y1 = 0; $y1 < $max; $y1++) { # go down the rows
            $y2 = $y1 + 1;
            $alt1 = ${$raxy}[$x1][$y1]; # LT
            $alt2 = ${$raxy}[$x2][$y1]; # RT
            $alt3 = ${$raxy}[$x2][$y2]; # RB
            $alt4 = ${$raxy}[$x1][$y2]; # LB
            $match = 0;
            $xm = $x1;
            $ym = $y1;
            # LT - RT
            $tx = $x1;
            $ty = $y1;
            if ($alt1 > $alt2) {
                $altl = $alt2;
                $alth = $alt1;
            } else {
                $altl = $alt1;
                $alth = $alt2;
            }
            if (($alt >= $altl) && ($alt <= $alth)) {
                ($lat,$lon) = get_latlon_xy($g_rmm_tile,$tx,$ty,$sqr);
                prt("$lon $lat ; $alt1 $alt2 ($tx,$ty)\n") if (VERB9());
                if ($match == 0) {
                    $xm = $tx;
                    $ym = $ty;
                }
                $match++;
            }
            # RT - RB
            $tx = $x2;
            $ty = $y1;
            if ($alt2 > $alt3) {
                $altl = $alt3;
                $alth = $alt2;
            } else {
                $altl = $alt2;
                $alth = $alt3;
            }
            if (($alt >= $altl) && ($alt <= $alth)) {
                ($lat,$lon) = get_latlon_xy($g_rmm_tile,$tx,$ty,$sqr);
                prt("$lon $lat ; $alt2 $alt3 ($tx,$ty)\n") if (VERB9());
                if ($match == 0) {
                    $xm = $tx;
                    $ym = $ty;
                }
                $match++;
            }
            # RB - LB
            $tx = $x2;
            $ty = $y2;
            if ($alt3 > $alt4) {
                $altl = $alt4;
                $alth = $alt3;
            } else {
                $altl = $alt3;
                $alth = $alt4;
            }
            if (($alt >= $altl) && ($alt <= $alth)) {
                ($lat,$lon) = get_latlon_xy($g_rmm_tile,$tx,$ty,$sqr);
                prt("$lon $lat ; $alt3 $alt4 ($tx,$ty)\n") if (VERB9());
                if ($match == 0) {
                    $xm = $tx;
                    $ym = $ty;
                }
                $match++;
            }
            # LB - LT
            $tx = $x1;
            $ty = $y2;
            if ($alt4 > $alt1) {
                $altl = $alt1;
                $alth = $alt4;
            } else {
                $altl = $alt1;
                $alth = $alt4;
            }
            if (($alt >= $altl) && ($alt <= $alth)) {
                ($lat,$lon) = get_latlon_xy($g_rmm_tile,$tx,$ty,$sqr);
                prt("$lon $lat ; $alt3 $alt4 ($tx,$ty)\n") if (VERB9());
                if ($match == 0) {
                    $xm = $tx;
                    $ym = $ty;
                }
                $match++;
            }
            # first match
            if ($match >= 2) {
                $tx = $xm;
                $ty = $ym;
                ($lat,$lon) = get_latlon_xy($g_rmm_tile,$tx,$ty,$sqr);
                prt("$lon $lat ; $alt1 $alt2 $alt3 $alt4 ($tx,$ty) $match\n"); #if (VERB9());
                ##$xg .= "$lon $lat ; $alt1 $alt2 $alt3 $alt4 ($tx,$ty) $match\n";
                ##$xg .= "NEXT\n";
                push(@matches,[$tx,$ty]);
                #$fnd = 1;
                #last;
            }
        }
        last if ($fnd);
    }
    my ($ra,$i,$j);
    $max = scalar @matches;
    for ($i = 0; $i < $max; $i++) {
        $ra = $matches[$i];
        $x1 = ${$ra}[0];
        $y1 = ${$ra}[1];
        $alt1 = ${$raxy}[$x1][$y1];
        ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x1,$y1,$sqr);
        $xg .= "$lon $lat ; $alt1\n";
        # find closest point
        for ($j = 0; $j < $max; $j++) {
            next if ($i == $j);
            $ra = $matches[$j];
            $x2 = ${$ra}[0];
            $y2 = ${$ra}[1];
            if (($x1 == $x2) && ($y1 == $y2)) {
                next;   # same point
            }
            if ($x1 == $x2) {
                if (abs($y1 - $y2) == 1) {
                    $alt2 = ${$raxy}[$x2][$y2];
                    ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x2,$y2,$sqr);
                    $xg .= "$lon $lat ; $alt2\n";
                }
            }
            if ($y1 == $y2) {
                if (abs($x1 - $x2) == 1) {
                    ($lat,$lon) = get_latlon_xy($g_rmm_tile,$x1,$y1,$sqr);
                    $alt2 = ${$raxy}[$x2][$y2];
                    $xg .= "$lon $lat ; $alt2\n";
                }
            }
        }
        $xg .= "NEXT\n";

    }
    return $max;
}


#########################################
### MAIN ###
parse_args(@ARGV);
##get_slippy_tile($def_lat,$def_lon,$def_zoom);
# use a 250x250 elevations file 
#process_in_file($def_file2);
#process_in_file($def_file3);
process_in_file($def_file4);

#process_in_file($in_file);
write_xg();
pgm_exit(0,"");
########################################

sub need_arg {
    my ($arg,@av) = @_;
    pgm_exit(1,"ERROR: [$arg] must have a following argument!\n") if (!@av);
}

sub parse_args {
    my (@av) = @_;
    my ($arg,$sarg);
    my $verb = VERB2();
    while (@av) {
        $arg = $av[0];
        if ($arg =~ /^-/) {
            $sarg = substr($arg,1);
            $sarg = substr($sarg,1) while ($sarg =~ /^-/);
            if (($sarg =~ /^h/i)||($sarg eq '?')) {
                give_help();
                pgm_exit(0,"Help exit(0)");
            } elsif ($sarg =~ /^v/) {
                if ($sarg =~ /^v.*(\d+)$/) {
                    $verbosity = $1;
                } else {
                    while ($sarg =~ /^v/) {
                        $verbosity++;
                        $sarg = substr($sarg,1);
                    }
                }
                $verb = VERB2();
                prt("Verbosity = $verbosity\n") if ($verb);
            } elsif ($sarg =~ /^l/) {
                if ($sarg =~ /^ll/) {
                    $load_log = 2;
                } else {
                    $load_log = 1;
                }
                prt("Set to load log at end. ($load_log)\n") if ($verb);
            } elsif ($sarg =~ /^o/) {
                need_arg(@av);
                shift @av;
                $sarg = $av[0];
                $out_file = $sarg;
                prt("Set out file to [$out_file].\n") if ($verb);
            } else {
                pgm_exit(1,"ERROR: Invalid argument [$arg]! Try -?\n");
            }
        } else {
            $in_file = $arg;
            prt("Set input to [$in_file]\n") if ($verb);
        }
        shift @av;
    }

    if ($debug_on) {
        prtw("WARNING: DEBUG is ON!\n");
        if (length($in_file) ==  0) {
            $in_file = $def_file;
            prt("Set DEFAULT input to [$in_file]\n");
        }
    }
    if (length($in_file) ==  0) {
        pgm_exit(1,"ERROR: No input files found in command!\n");
    }
    if (! -f $in_file) {
        pgm_exit(1,"ERROR: Unable to find in file [$in_file]! Check name, location...\n");
    }
}

sub give_help {
    prt("$pgmname: version $VERS\n");
    prt("Usage: $pgmname [options] in-file\n");
    prt("Options:\n");
    prt(" --help  (-h or -?) = This help, and exit 0.\n");
    prt(" --verb[n]     (-v) = Bump [or set] verbosity. def=$verbosity\n");
    prt(" --load        (-l) = Load LOG at end. ($outfile)\n");
    prt(" --out <file>  (-o) = Write output to this file.\n");
}

# eof - template.pl

index -|- top

checked by tidy  Valid HTML 4.01 Transitional