#!/usr/bin/perl
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.html
# Created: 2024-02-05T23:40:23+0000
# Last-Updated: 2024-05-16T01:01:22+0000
#     Update #: 207

=head1 NAME

PLSS_projection - Using projections on PLSS data

=head1 DESCRIPTION

We place an oblique Mercator projection on each individual square mile
PLSS section. allowing queries, forward and inverse.

We should use the southEAST corner as our origin, as this is usually the
most reliable (north tops of sections are cut off more often, and west
edges of sections adsorb surveying errors often at township edges.)

However, as all these programs I am writing are just demonstration
programs, we shall instead take advantage of the PID _xxxyyy format, which
uses the southWEST corner of each section as the origin. Thus at the
gain of much simpler "what is the id of the origin point for this
section" calculations, we lose some occasional precision.

We lay our "equator" upon the south edge of the section, and get the
azimuth of the SW corner to the SE corner.

Now with ease we can compute where e. g., "one city block north and
east of the origin section corner" is, by asking for _x10y10 .

=cut

package PLSS_projection;
use strict;
use warnings q(all);
use Bou2colrow;
use POSIX;
use Math::Trig qw(:great_circle deg2rad rad2deg);
use Data::Dumper;
sub NESW { deg2rad( $_[0] ), deg2rad( 90 - $_[1] ) }
$Data::Dumper::Indent   = 1;
$Data::Dumper::Sortkeys = 1;

=head2 pid2lonlat

Given a corners file,

 -87.8677047,42.1529536,IL030420N0120E0_200700
 -87.8680876,42.1383760,IL030420N0120E0_200600
 -87.8874919,42.1384733,IL030420N0120E0_100600...

And a targets file,

 IL030420N0120E0_210610,,3500 (W etc.)
 (the empty field is for compatibility with a pid2addr
 program...)

We then output each target's coordinates.

 -87.865661,42.140192,3500
 -87.863227,42.140197,3400

=cut

sub pid2lonlat {
    die "Usage: need two files, like corners.csv targets.csv" unless @ARGV == 2;
    my %twsp;
    open( my $corners, "<", $ARGV[0] )
      or die "Can't open corners file $ARGV[0]: $!";
    while (<$corners>) {
        chomp;
        my @F   = split /,/;
        my @pid = split /_/, $F[2];
        $twsp{ $pid[0] }{ $pid[1] }{coordinates} = [ @F[ 0, 1 ] ]
          ;    #absorbs duplicates too, e.g., adjacent sections' shared corners.
    }
    close $corners;

    ## Compute azimuths
    for my $t ( keys %twsp ) {
        for my $p ( keys %{ $twsp{$t} } ) {
            $p =~ /^([1-7])00([1-7])00$/ or die "Odd $_";
            next unless $1 < 7;
            my $east = $1 + 1 . 0 . 0 . $2 . 0 . 0;
            next unless $twsp{$t}{$east};
            my @L = NESW( @{ $twsp{$t}{$p}{coordinates} } );
            my @T = NESW( @{ $twsp{$t}{$east}{coordinates} } );
            $twsp{$t}{$p}{azimuth} = rad2deg great_circle_direction( @L, @T );
        }
    }
    open( my $targets, "<", $ARGV[1] )
      or die "Can't open targets file $ARGV[1]: $!";
    my %section_status;
    my @output;
    my @re = (
        qr/^[A-Z]{2}\d{2}\d{3}0[NS]\d{3}0[EW]0$/,
        qr/^[1-7][0-7][0-9][1-7][0-7][0-9]$/
    );
    while (<$targets>) {
        chomp;
        my @F   = split /,/;
        my @pid = split /_/, $F[0];
        for ( 0, 1 ) {
            die "\"$pid[$_]\" !~ \"$re[$_]\"" unless $pid[$_] =~ $re[$_];
        }
        my @xxxyyy  = split //, $pid[1];
        my $SWid    = $xxxyyy[0] . 0 . 0 . $xxxyyy[3] . 0 . 0;
        my $azimuth = $twsp{ $pid[0] }{$SWid}{azimuth};
        my $PIDSNSW = Bou2colrow::xxxyyy2SN_US( $F[0] ) . "[$SWid]";
        unless ( defined $azimuth ) {
            $section_status{$PIDSNSW}--;
            $section_status{badlines}{$_}++;
            next;
        }
        my @miles;
        for ( 0, 1 ) {
            my $offset = 3 * $_;
            push @miles,
              $xxxyyy[ 1 + $offset ] / 8 + $xxxyyy[ 2 + $offset ] / 80;
        }
        my $cmd =
          sprintf
          "echo @miles $F[2] | proj -f %%.6f -I +proj=omerc +alpha=%s +gamma=0 "
          . "+lonc=%s +lat_0=%s +units=us-mi\n",
          90 - $twsp{ $pid[0] }{$SWid}{azimuth},
          @{ $twsp{ $pid[0] }{$SWid}{coordinates} };
        my $o = qx/$cmd/;  #one day use IPC::Run3 when it becomes more standard.
        for ($o) { s/ /\t/; s/\t/,/g; print; }
        $section_status{$PIDSNSW}++;
    }
    close $targets;
    print STDERR
      "Sections [and SW corner ids] use count summary, just for your reference.$/";
    print STDERR Data::Dumper->Dump( [ \%section_status ],
        [qw(*section_status)] );
    die "Get all missing first, please. They are the negative numbers above."
      if $section_status{badlines};
}

1;