Hello, all--
I had a friend who wanted to find some survey markers near his house, so I looked around and found a webservice that would locate some survey markers in an area. Since I coded it up, I thought I'd publish it here, in case there's anyone else who might like to try it. It's rough (as it's a one-off), but should be easy enough to modify.
As is often the case, all the heavy lifting is done by some handy CPAN modules (JSON, HTTP::Request, LWP::UserAgent and Math::Trig). I was especially pleased to find Math::Trig--I was trying to derive it myself and needed arc-cosine. When I found Math::Trig had acos and looked over the docs, I found that it already had all the great-circle math as well!
Anyway, I hope someone finds it useful...
#!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 lati
+tude
# 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 ci
+rcle)
EOMSG
my $Re_mi = 3958.8; # radius of earth in miles
# Figure how approximately how long a degree is in both the longitudin
+al 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 sea
+rch
# 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'=>'applica
+tion/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 cen
+tered 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%02
+d.%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 <x
+DDDMMSS.s*>
if ($dir eq "N" or $dir eq "S") {
return sprintf "%s%02d%02d%02d.%s", $dir, $deg, $min, $sec, $s
+frac;
}
else {
return sprintf "%s%03d%02d%02d.%s", $dir, $deg, $min, $sec, $s
+frac;
}
}
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 ];
}
A sample run gives me:
$ perl websvc_usgs_fetch_bounding_box.pl 34.1234 -78.1234 1.5
Found 11 markers within 1.5 miles centered around [-78.1234, 34.1234]
LATITUDE LONGITUDE NAME PID
---------- ---------- -------------------- ----------
34.14829 -78.09827 WINNABOW RM 4 EB0192
34.14832 -78.09805 WINNABOW EB0189
34.14856 -78.09788 WINNABOW RM 3 EB0191
34.09944 -78.12361 P 117 EB0198
34.11095 -78.11490 FLOWERS EB2124
34.11178 -78.11387 Q 117 EB0197
34.11812 -78.10776 BECK EB2125
34.12278 -78.10611 R 117 EB0195
34.12708 -78.10360 CAMPBELL EB2108
34.14593 -78.10423 WINNABOW AZ MK 2 EB1389
34.14722 -78.11694 N 235 EB0300
...roboticus
When your only tool is a hammer, all problems look like your thumb.
-
Are you posting in the right place? Check out Where do I post X? to know for sure.
-
Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:
<code> <a> <b> <big>
<blockquote> <br /> <dd>
<dl> <dt> <em> <font>
<h1> <h2> <h3> <h4>
<h5> <h6> <hr /> <i>
<li> <nbsp> <ol> <p>
<small> <strike> <strong>
<sub> <sup> <table>
<td> <th> <tr> <tt>
<u> <ul>
-
Snippets of code should be wrapped in
<code> tags not
<pre> tags. In fact, <pre>
tags should generally be avoided. If they must
be used, extreme care should be
taken to ensure that their contents do not
have long lines (<70 chars), in order to prevent
horizontal scrolling (and possible janitor
intervention).
-
Want more info? How to link
or How to display code and escape characters
are good places to start.