#!/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 = ; 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 = ; 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 = < $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 (-o) = Write output to this file.\n"); } # eof - template.pl