#!/usr/bin/perl
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.htm
# Created: 2023-11-20T22:22:01+0000
# Last-Updated: 2024-05-16T02:34:37+0000
#     Update #: 457
#

=head1 NAME

Twsp2mi - Convert PLSS Townships/Ranges+sections point ids to miles and back

=head1 DESCRIPTION

Here we dissect "Bureau of Land Management Geographic Coordinate
Database (GCDB) Publication Point IDs", e.g., AZ230490N0130W0_240400 .

See also "GCDB NAMING CONVENTIONS Meridian codes, Township names,
Point IDs, SID Numbers".

Note we expect everything is within the same meridian (but interstate is OK.)

=cut

use strict;
use warnings q(all);

package Twsp2mi;
use Bou2colrow;

=head1 mi_diff

Given pointId0, pointId1, return diff_in_milesEW, diff_in_milesNS;

 my @m = mi_diff(
     qw!IL030010N0010E0_100100
        IL030430N0090E0_340105!
 );
 print "@m\n";

 prints
 50.5 252.0625

Safe crossing to different sides of meridians and baselines,

 my @m = mi_diff(
     qw!
       IL030010N0010E0_100100
       IL030010S0090E0_340105!
 );
 print "@m\n";
 #50.5 -5.9375

=cut

sub mi_diff {

    my ( @i, @miles );
    die "Need two arguments" unless $#_ == 1;
    for (@_) {
        push @i, [
            unpack
              ## AZ23 049 0 N 013 0 W 0 _ 2 4 0 4 0 0
              "  A2A2  A3 A A  A3 A A A A A A A A A A",
            ##    0 1   2 3 4   5 6 7 8 9 0 1 2 3 4 5
            $_,
        ];
    }
    die "Not ready for different meridians"
      unless $i[0][1] == $i[1][1];

    ## We didn't use our state check parsing... and in fact we shouldn't,
    ## as some meridians are valid interstate.

    ## could also check for TxW TxE RxN RxS but maybe Ohio...
    die "No '_' in both @_"
      unless $i[0][9] eq $i[1][9] eq "_";
    for my $field ( 3, 6, 8 ) {
        die "Not ready for non-zero at item $field + 1 in
        @_"
          unless $i[0][$field] eq $i[1][$field] eq 0;
    }
    my @L =
      ( [ 5, 7, 10, 11, 12 ], [ 2, 4, 13, 14, 15 ] );
    for ( 0, 1 ) {
        push @miles,
          twdiffi(
            [ @{ $i[0] }[ @{ $L[$_] } ] ],
            [ @{ $i[1] }[ @{ $L[$_] } ] ],
          );
    }
    return @miles;
}

=head1 twdiffi

(Retired for twdiffi_chains. See below.)

Computes the difference in miles, in one dimension.

 print twdiffi( [ 165, 'W', 1, 0, 0 ], [ 164, 'W', 1, 4, 5 ] );    #6.5625
 # Fields are [$townships, $nesw, $rows, $furlongs, $chains].

Yes, can cross from one side of a meridian or baseline to the other safely,

 print twdiffi( [ qw/1 W 6 4 5/ ], [ qw/1 E 1 0 0/] );
 #prints 0.4375

=cut

sub Otwdiffi {
    my $total;
    for ( 0, 1 ) {
        my ( $townships, $nesw, $rows, $furlongs, $chains ) = @{ $_[$_] };
        ## T1S100 and R1W100 are six miles south or west of the meridian or baseline:
        if ( $nesw =~ /S|W/ ) { $townships *= -1; }
        ## But T1N100 and R1E100 are right on the meridian or baseline:
        else { $townships -= 1 }
        my $miles = $townships * 6;
        $miles  += $rows - 1;
        $chains += $furlongs * 10;
        $miles +=
          $chains / 80;    # A jidanni extension. Before only 0 and 5 were OK.
        $total += $miles * ( $_ ? 1 : -1 );    #subtract FROM from TO.
    }
    return $total;
}

=head1 twdiffi_chains

A safer version of twdiffi. Avoids all fractions.

=cut

sub twdiffi_chains {
    my $total;
    for ( 0, 1 ) {
        my ( $townships, $nesw, $rows, $furlongs, $chains ) = @{ $_[$_] };
        ## T1S100 and R1W100 are six miles south or west of the meridian or baseline:
        if ( $nesw =~ /S|W/ ) { $townships *= -1; }
        ## But T1N100 and R1E100 are right on the meridian or baseline:
        else { $townships -= 1 }
        my $miles = $townships * 6;
        $miles  += $rows - 1;
        $chains += $furlongs * 10;
        $chains += $miles * 80;
        $total  += $chains * ( $_ ? 1 : -1 );    #subtract FROM from TO.
    }
    return $total;
}

=head1 point_id2miles

(Obsolete: use point_id2chains)

Given an id IL03..._... etc., return both how many miles (east is
positive) away from the (NS) meridian, and away (north is positive)
from the (EW) baseline.

 for (qw!IL030010N0010E0_100100 IL030430N0090E0_340105!) {
     my @m = point_id2miles($_); #Obsolete. Use point_id2chains.
     print "$_: @m\n";
 }

 prints:
 IL030010N0010E0_100100: 0 0
 IL030430N0090E0_340105: 50.5 252.0625

=cut

sub Opoint_id2miles {
    my @ans;
    for ( 0, 1 ) {
        my ( $townships, $nesw, $rows, $furlongs, $chains ) = pid2extract( @_, $_ );
        ## Gosh, the following reuses code from above. I really should make a third sub{}.
        ## T1S100 and R1W100 are six miles south or west of the meridian or baseline:
        if ( $nesw =~ /S|W/ ) { $townships *= -1; }
        ## But T1N100 and R1E100 are right on the meridian or baseline:
        else { $townships -= 1 }
        my $miles = $townships * 6;
        $miles  += $rows - 1;
        $chains += $furlongs * 10;
        $miles +=
          $chains / 80;  # A jidanni extension. Before only 0 and 5 were OK.
                         #The above / is the source of floating point artifacts.
            #I should do all math in eightieths of a mile whole numbers.
        push @ans, $miles;
    }
    return @ans;
}

=head1 point_id2chains

A safer version of point_id2miles, avoids all fractions.

=cut

sub point_id2chains {
    my @ans;
    for ( 0, 1 ) {
        my ( $townships, $nesw, $rows, $furlongs, $chains ) = pid2extract( @_, $_ );
        ## Gosh, the following reuses code from above. I really should make a third sub{}.
        ## T1S100 and R1W100 are six miles south or west of the meridian or baseline:
        if ( $nesw =~ /S|W/ ) { $townships *= -1; }
        ## But T1N100 and R1E100 are right on the meridian or baseline:
        else { $townships -= 1 }
        my $miles = $townships * 6;
        $miles  += $rows - 1;
        $chains += $furlongs * 10;
        ## Before, $chains could only 0 and 5. But I, jidanni hereby
        ## have extended it, ignoring other uses of that digit.
        $chains += $miles * 80;
        push @ans, $chains;
    }
    return @ans;
}

=head1 pid2extract

Given e.g., IL03..._... , get TR, NESW, rows, furlongs, chains. Yes,
furlongs are eighths of a mile.

Originally Chains can only be 0, or 5 (half a furlong) but I'll allow
Chains to be just are eightieths of a mile. My non-standard
enhancement...

Second argument: 0 for EW, 1 for NS.

 for ( 0, 1 ) {
     my @m = pid2extract( "IL030430N0090E0_340105", $_ );
     print "@m\n";
 }

 prints
 009 E 3 4 0
 043 N 1 0 5

Note we don't check what's in non-standard fields distinguishing
identical townships, etc.

=cut

sub pid2extract {
    die "No defined argument passed to me"           unless defined $_[0];
    die "Need two aruments"           unless @_ == 2;
    die "Second argument must be [01]" unless $_[1] =~ /^[01]$/;
    die "\"$_[0]\" is not a standard ..._... length" unless length $_[0] == 22;
    ## Garbage input shorter than needed will cause perl to
    ## "raise an exception."
    my @K = unpack $_[1]    #0: R:EW; 1: T:NS
      ##   AZ23 049 0 N 013 0 W 0 _ 2 4 0 4 0 0
      ? "xxxx  A3 x A xxx x x x x x x x A A A"
      : "xxxx xxx x x  A3 x A x x A A A x x x",
      $_[0];
    my $field;
    $field = 2;
    for ( $K[$field] ) {
        if ( $_ < 1 || $_ > 7 ) {
            er( \@K, $field );
        }
    }
    $field = 3;
    for ( $K[$field] ) {
        if ( $_ > 7 ) {
            er( \@K, $field );
        }
    }
    return @K;

    sub er {
        die sprintf "@{$_[0]}: invalid field %d/%d: $_", $_[1],
          scalar @{ $_[0] };
    }

}

=head1 snid2extract()

Like pid2extract but for ...SN... ids. 2nd argument is NW, NE, C etc. corner.

=cut

sub snid2extract{
    return pid2extract(Bou2colrow::trsec_dir2id(@_));
}

=head1 Functions to create IDs from miles / chains:

(Chains (mile/80) needed to avoid fractional arithmetic residues.)

=head1 chains2id1d

Here given chains from baseline or meridian, compute the components,
either the X or the Y of QQ00YYY0YXXX0X0_XXXYYY. And return them as an
array e.g., (YYY Y YYY), where MN460280N0220W0_411322 Y's would be (28
1 322)=(T or R 28, N or E, _ 322yyy or xxx322). 1 is N or E, -1 is S
or W.

 for ( -3 .. 3 ) { my @o = chains2id1d($_); print "@o$/"; }
 #prints:
 1 -1 6 7 7
 1 -1 6 7 8
 1 -1 6 7 9
 1 1 1 0 0
 1 1 1 0 1
 1 1 1 0 2
 1 1 1 0 3

 for ( 23, -23, -8 * 80 - 1 ) {
     my @z = chains2id1d($_);
     print qq(@z), $/;
 }
 #prints:
 1 1 1 2 3
 1 -1 6 5 7
 2 -1 4 7 9

=cut

sub chains2id1d {
    my $chains    = $_[0];
    my $sign      = $chains < 0 ? -1 : 1;
    my $townships = int abs( $chains / ( 80 * 6 ) ) + 1;    #No T or R 0!
    my $township_chains =
      ( $chains + 80 * 99999 * 6 ) % ( 80 * 6 )
      ;    # chains within one township, counted always from W or S edges.
    my $sections       = int( $township_chains / 80 ) + 1;       #1..6, not 0..5
    my $furlongs       = int( ( $township_chains % 80 ) / 10 );  #0..7
    my $furlong_chains = $township_chains %
      10;    #chains within one furlong, 0..9 (was only 0, 5 but I extended it!)
    return $townships, $sign, $sections, $furlongs, $furlong_chains;
}

=head1 point_id_plus_chains

Given a point id, and CHAINS +-east, and +-north of it to add, return a point id with them added

(We need to use CHAINS (1/80 of a mile,) to avoid fractions. Giving
this function a non-integer of chains will trip an alarm. (Dies with error message.)

 print point_id_plus_chains( 'MN460280N0220W0_440300', -1, -2 ), $/;
 #prints MN460280S0220E0_439278

=cut

sub point_id_plus_chains {
    die "Wrong number of arguments in: @_" unless @_ == 3;
    my @chains = point_id2chains( $_[0] );
    my @r;
    for ( 0, 1 ) {
        for ( $_ + 1 ) {
            die "$_[$_] is not an integer number of chains!"
              if $_[$_] - int $_[$_];
        }
        $chains[$_] += $_[ $_ + 1 ];
        push @r, [ chains2id1d( $chains[$_] ) ];
    }

    return sprintf "%s%03d0%s%03d0%s0_%d%d%d%d%d%d",
      substr( $_[0], 0, 4 ),
      $r[1][0],
      $r[1][1] < 0 ? "S" : "N",
      $r[0][0],
      $r[0][1] < 0 ? "W" : "E",
      @{ $r[0] }[ 2, 3, 4 ],
      @{ $r[1] }[ 2, 3, 4 ];
}

=head1 point_id_plus_miles

Wrapper for the above.

 print point_id_plus_miles( 'MN460280N0220W0_440300', -1.5, -2.25 ), $/;
 MN460270N0220W0_300660

Yes, point_id_plus_miles( 'MN460280N0220W0_440300', -1.5, -2.2587 ) will trigger our
non-integer number of chains burglar alarm systems that live in point_id_plus_chains():
"-180.696 is not an integer number of chains!"

=cut

sub point_id_plus_miles {
    return point_id_plus_chains( $_[0], $_[1] * 80, $_[2] * 80 );
}

=head1 pid2gcp_mi

Given lines with PIDs, return -gcp lines to be used on the
gdaltransform [--optfile] command line, now having miles from PLSS
origin.

 $ cat gcp_seed.csv
 -87.8874918582661,42.1384732962006,IL030420N0120E0_100600
 -87.8880983169846,42.0804153949808,IL030420N0120E0_100100
 -87.7712238884793,42.0797003951774,IL030420N0120E0_700100
 $ perl -MTwsp2mi -we Twsp2mi::pid2gcp_mi gcp_seed.csv
 -gcp 66 251 -87.8874918582661 42.1384732962006
 -gcp 66 246 -87.8880983169846 42.0804153949808
 -gcp 72 246 -87.7712238884793 42.0797003951774

=cut

sub pid2gcp_mi {
    while (<>) {
        chomp;
        my @F = split /,/;
        ## Obsolete. Need to change to point_id2chains:
        print join " ", "-gcp", point_id2miles( $F[2] ), @F[ 0, 1 ];
        print $/;
    }
}

1;