#!env perl
#
# websvc_usgs_fetch_bounding_box.pl  <LAT> <LON> <dist>
#
# Use the USGS "Bounding Box Service" to find survey markes within the
# rough rectangle whose sides are <dist> miles from the specified latitude
# and longitude.
#
# 20180504 original version
#
use strict;
use warnings;
use Data::Dump 'pp';
use HTTP::Request;
use JSON;
use LWP::UserAgent;
use Math::Trig qw( :great_circle deg2rad );

my $LAT = shift;
my $LON = shift;
my $center_dec = [ $LON, $LAT ];
my $squaradius_mi = shift or die <<EOMSG;
Expected:
perl websvc_usgs_fetch_bounding_box.pl LAT LON RAD

    LAT - latitude like 38.1234,
    LON - longitude like -78.1234,
    RAD - radius in miles (actually roughly a rectangle rather than circle)
EOMSG

my $Re_mi = 3958.8; # radius of earth in miles

# Figure how approximately how long a degree is in both the longitudinal and
# latitudinal directions.
my $mi_per_degree = miles_per_degree([ NESW($LON, $LAT) ]);

# Now find the (min/max) * (lat,lon) for the bounding rectangle to search
# for survey markers
my $deg_per_mi_lat = 1 / $mi_per_degree->[1];
my $deg_per_mi_lon = 1 / $mi_per_degree->[0];

my $min_lat = $center_dec->[1] - $squaradius_mi*$deg_per_mi_lat;
my $max_lat = $center_dec->[1] + $squaradius_mi*$deg_per_mi_lat;
my $min_lon = $center_dec->[0] - $squaradius_mi*$deg_per_mi_lat;
my $max_lon = $center_dec->[0] + $squaradius_mi*$deg_per_mi_lat;

# We'll use the DDMMSS format for the lat/lon
$min_lat = dec_to_DDMMSS($min_lat, "N", "S");
$max_lat = dec_to_DDMMSS($max_lat, "N", "S");
$min_lon = dec_to_DDMMSS($min_lon, "E", "W");
$max_lon = dec_to_DDMMSS($max_lon, "E", "W");

# Now fetch the data, and print the results
my $URL = "http://geodesy.noaa.gov/api/nde/bounds?"
    ."minlat=$min_lat&maxlat=$max_lat"
    ."&minlon=$min_lon&maxlon=$max_lon";

my $request = HTTP::Request->new(GET=>$URL, [ 'Content-Type'=>'application/json; charset=UTF-8' ]);
my $ua = LWP::UserAgent->new;
my $response = $ua->request($request);

if (! exists $response->{_content}) {
    # crappy error detection/handling but it meets my current needs
    print pp($response), "\n";
    print "***** expected response->{content}!!!!!!!!\n";
}

my $r = decode_json($response->{_content});

print "Found ", scalar(@$r), " markers within $squaradius_mi miles centered around ",
      pp($center_dec), "\n";

print <<EOHDR;
LATITUDE   LONGITUDE  NAME                 PID
---------- ---------- -------------------- ----------
EOHDR
for my $hr (@$r) {
    printf "%-10s %-10s %-20s %-10s\n", $hr->{lat}, $hr->{lon}, $hr->{name}, $hr->{pid};
}

sub dec_to_DDMMSS {
    my ($dec, $dir_pos, $dir_neg) = @_;

    # Unfortunately, they use 2 sig digs for N/S and 3 for E/W
    my $fmt = $dir_pos eq 'N' ? "%s%02d%02d%02d.%03d" : "%s%03d%02d%02d.%03d";

    my $dir = $dir_pos;
    if ($dec < 0) { $dec = -$dec; $dir = $dir_neg; }

    my ($deg, $min, $sec, $sfrac) = (int($dec), 0, 0, 0);

    $dec = 60 * ($dec - $deg);
    $min = int($dec);
    $dec = 60 * ($dec - $min);
    $sec = int($dec);
    $sfrac = int(1000*($dec-$sec));

    return sprintf $fmt, $dir, $deg, $min, $sec, $sfrac;
}

sub dms_to_DDMMSS {
    my $dms = shift;
    my ($deg, $dir, $min, $sec, $sfrac);
    if ($dms =~ /^\s*(\d+)\s*([NEWS])\s*(\d+)\s*'\s*([\d\.]+)\s*"\s*$/) {
        ($deg, $dir, $min, $sec, $sfrac) = ($1, $2, $3, $4, 0);
        if ($sec =~ /(\d+)\.(\d+)/) {
            ($sec,$sfrac) = ($1,$2);
        }
    }
    else {
        die "Unexpected format <$dms>!";
    }

    # Build the return value: For N/S use <xDDMMSS.s*>, for E/W use <xDDDMMSS.s*>
    if ($dir eq "N" or $dir eq "S") {
        return sprintf "%s%02d%02d%02d.%s", $dir, $deg, $min, $sec, $sfrac;
    }
    else {
        return sprintf "%s%03d%02d%02d.%s", $dir, $deg, $min, $sec, $sfrac;
    }
}

sub NESW { deg2rad($_[0]), deg2rad(90-$_[1]) }

sub dms_to_dec {
    my $dms = shift;
    my ($deg, $dir, $min, $sec, $sfrac);

    if ($dms =~ /^\s*(\d+)\s*([NEWS])\s*(\d+)\s*'\s*([\d\.]+)\s*"\s*$/) {
        ($deg, $dir, $min, $sec, $sfrac) = ($1, $2, $3, $4, 0);
        if ($sec =~ /(\d+)\.(\d+)/) {
            ($sec,$sfrac) = ($1,$2);
        }
    }
    else {
        die "Unexpected format <$dms>!";
    }

    $dir = ($dir eq 'N' or $dir eq 'E') ? 1 : -1;
    return $dir * ($deg + $min/60.0 + (0 + ("$sec.$sfrac"))/3600.0);
}

# compute number of miles per degree at the specified lat/lon
sub miles_per_degree {
    my $news = shift;

    my $news_1lat = [ @$news ];
    $news_1lat->[1] += deg2rad(1.0);

    my $news_1lon = [ @$news ];
    $news_1lon->[0] += deg2rad(1.0);

    my $dLat_km = great_circle_distance(@$news, @$news_1lat, $Re_mi);
    my $dLon_km = great_circle_distance(@$news, @$news_1lon, $Re_mi);

    return [ $dLat_km, $dLon_km ];
}