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.
Back to
Cool Uses for Perl