#!/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-04-27T00:34:03+0000
#     Update #: 558

=head1 NAME

PointId2Address - Turn US PLSS point ids into address values, or vice versa, based on parameters given

=cut

use strict;
use warnings q(all);
use Data::Dumper;

package PointId2Address;
## use Data::Dump qw(dump);
## $Data::Dump::LINEWIDTH = 11;

use Twsp2mi;

=head1 SYNOPSIS

Functions:

=head1 twsp_blanket

 print "$_\n" for
  twsp_blanket(
    {
        state => "IL",
        meridian => "03",
        T    => [ -35 .. -26 ],
        R    => [ -65 .. -52 ],
        S    => [ 31, 33, 35 ]
    }
  );

returns all their combinations, as SN plssids.

 IL030350S0650W0SN310
 IL030350S0650W0SN330
 IL030350S0650W0SN350
 ...

(Negatives: -42 = T42S or R42W).

Can of course also be used to get a strip, one township wide, e.g., by selecting only one Range.

Or just pretty print a single township's single section,

 print "$_\n"
   for twsp_blanket(
     {
         state    => "IL",
         meridian => "03",
         T        => [40],
         R        => [13],
         S        => [23],
     }
   );

gives IL030400N0130E0SN230

We can also produce just townships, no sections:

 print "$_\n"
   for PointId2Address::twsp_blanket(
     {
	 state    => "SD",
	 meridian => "05",
	 T        => [128],
	 R        => [ -80 .. -47 ],
     }
   );

 SD051280N0800W0
 SD051280N0790W0
 SD051280N0780W0...

And we can get intersecting strips. Hard to explain...

 my %h = (
     t => [ [23], [ -12 .. 23 ] ],    #T23N and T12S to T23N
     r => [ [1],  [ 1 .. 29 ] ]       #R1E and R1E to R29E
 );
 my %check
   ; #A hash, so we can remove the doubled node at the origin where the two axes meet.
 my @wants;    #But we like our array order vs. hash's...
 for ( 0, 1 ) {
     for (
	 PointId2Address::twsp_blanket(
	     {
		 state    => "SD",
		 meridian => "07",           #Black Hills
		 T        => $h{t}[$_],
		 R        => $h{r}[ !$_ ],
	     }
	 )
     ) {
	 push @wants, $_ unless $check{$_}++;
     }
 }
 print "$_\n" for @wants;

 SD070230N0010E0
 SD070230N0020E0
 SD070230N0030E0...

=cut

sub twsp_blanket {
    use Data::Dump q!pp!;
    my @res;
    for my $T ( @{ $_[0]->{T} } ) {
        next unless $T;    #skip over "T0"
        my $t = sprintf "%03d0%s", abs $T, $T < 0 ? "S" : "N";
        for my $R ( @{ $_[0]->{R} } ) {
            next unless $R;    #skip over "R0"
            my $r   = sprintf "%03d0%s", abs $R, $R < 0 ? "W" : "E";
            my $str = sprintf "%s%s%s%s0",
              $_[0]->{state},
              $_[0]->{meridian},
              $t, $r;
            unless ( $_[0]->{S} ) {
                push @res, $str;
                next;
            }
            for my $s ( @{ $_[0]->{S} } ) {
                push @res, sprintf "%sSN%02d0", $str, $s;
            }
        }
    }
    return @res;
}

=head1 point_id_diff

Let's take Chicago (North side) for example.

How many miles up and over, is Howard and Western, from Addison and Pulaski?

Let's see. Addison & Pulaski is the half mile point between sections
22 and 23, T40N R13E IL03 meridian.

If we can't or don't want to figure it out on paper, we transform
IL030400N0130E0SN230 (that we just obtained above in the twsp_blanket
example,) via

 use Bou2colrow;
 print Bou2colrow::trsec_dir2id( qw!IL030400N0130E0SN230 W! );
 #getting IL030400N0130E0_500340

Also we use IL030410N0140E0_100240 for Howard & Western (obtained in example further down.)

 my @a = point_id_diff(
     {
         origin => {    #Addison & Pulaski
             id => "IL030400N0130E0_500340",
             ## miles_from_pm_bl => [ 76, 236.5 ],
         },
         target => {    #Howard & Western
             id => "IL030410N0140E0_100240",
             ## miles_from_pm_bl => [ 78, 241.5 ],
         },
     }
 );

miles_from_pm_bl: miles from principal meridian and baseline of this
id. Can pre-compute and pass it in for big batch jobs. Or some of our
programs cache it, so the user might want to clear it if needed.

OK, @a has our answer: (2, 5). Two miles east and four miles north.

Wait, the above is OK for one shot calls.
Let's go back and see how we can use the {miles} recomputing saving feature.
(Highlights of ../us/nd/north_dakota_state/burkle:)

 my %g;
 $g{origin}{id} = "ND051460N0780W0_700100";
 while (<>) {
     chomp;
     my @F = split /,/;
     $g{target} = { id => $F[2] };    #we wipe out the previous {target}{miles}
     my @a = PointId2Address::point_id_diff( \%g );
     ...

OK, we have safely got our results in @a.

=cut

sub point_id_diff {
    my %grid = %{ $_[0] };
    for (qw/target origin/) {
        next if $grid{$_}{miles_from_pm_bl};
        @{ $grid{$_}{miles_from_pm_bl} } =
          Twsp2mi::point_id2miles( $grid{$_}{id} );
    }
    for ( 0, 1 ) {
        $grid{diff}[$_] =
          $grid{target}{miles_from_pm_bl}[$_] -
          $grid{origin}{miles_from_pm_bl}[$_];
    }
    return @{ $grid{diff} };
}

=head1 id2addr

We start in Chicago at State and Madison (0, 0),

 IL030390N0140E0SN100 SW =
 IL030390N0140E0_400500

What might be the house numbers at Howard & Western? First let's get the point ID there,

 use Bou2colrow;
 Bou2colrow::trsec_dir2id( qw!IL030410N0130E0SN250 E! );
 #gives IL030410N0140E0_100240

 my @a = id2addr(
     {
         num_per_mile => [ 800, 800 ], # Always positive numbers
         origin       => {               #State & Madison
             address => [ 0, 0 ], #actually the default, if omitted.
             id      => "IL030390N0140E0_400500",
 ##          miles_from_pm_bl => [ 81, 232 ],
         },
         target => {                     #Howard & Western
             id => "IL030410N0140E0_100240",
 ##          miles_from_pm_bl => [ 78,  241.5 ],
         }
     }
 );

 # miles_from_pm_bl: see previous documentation.
 printf "%s%s %s%s\n", @a; #Gives 2400W 7600N, our answer!

How about Lake Co. Illinois, whose origin at the SW corner is not 0,0 but something else.
No problem, just plug it in...

 ...      origin => {
              address => [ -29000, 20000 ],...
          },
          num_per_mile => [ (1000) x 2 ] # Always positive numbers

=cut

sub id2addr {
    return EWNS2( id2addr_raw(@_) );
}

=head1 id2addr_raw()

Same as id2addr, but just returns the two numbers, e.g., -600, 800, no
N S E W jazz.

=cut

sub id2addr_raw {
    my %grid = %{ $_[0] };
    for (qw/target origin/) {
        die "Didn't get {$_} passed to me" unless $grid{$_};
        ## Maybe we can store this in a global variable in order to only compute once on big jobs.
        ## Or maybe perl does something behind the scenes to avoid recomputing the same number...
        next if $grid{$_}{miles_from_pm_bl};
        @{ $grid{$_}{miles_from_pm_bl} } =
          Twsp2mi::point_id2miles( $grid{$_}{id} );
    }
    for ( 0, 1 ) {
        $grid{diff}[$_] =
          $grid{target}{miles_from_pm_bl}[$_] -
          $grid{origin}{miles_from_pm_bl}[$_];
        $grid{target}{address}[$_] =
          ( $grid{target}{miles_from_pm_bl}[$_] -
              $grid{origin}{miles_from_pm_bl}[$_] ) * $grid{num_per_mile}[$_] +
          ( $grid{origin}{address}[$_] || 0 );
    }
    return @{ $grid{target}{address} };
}

=head1 addr2id

Given an address grid origin definition, and house numbers X an Y, e.g., -600, 100
returns e.g., MN460280N0220W0_415300.

 my $id = addr2id(
     {
	 num_per_mile => [ 1600, 800 ], # Always positive numbers
	 origin       => {
	     address => [ 0, 0 ],    #actually the default, if omitted.
	     id      => "MN460280N0220W0_440300",
	 },
     },
     -600,
     100
 );

 #result: $id = "MN460280N0220W0_415300"

Note many cities blocks go from 1200-1260, 1300-1360... with the half block at x30 not x50.

Also some cities don't have 000 blocks. I.e., 300-200S, 200-100S, 100-200N, 200-300N...

All that needs to be translated into "regular blocks" before being fed
in here, and also translated back out when returned. That is all the
responsibility of the calling functions.

=cut

sub addr2id {
    my %grid = %{ (shift) };
    my $info = qq(@_);
    my @a2;
    for ( 0, 1 ) {
        $a2[$_] = $_[$_] - @{ $grid{origin}{address} }[$_] || 0;
    }

    ## Cannot make miles, must make chains, to avoid fractions.
    my @chains;
    for ( 0, 1 ) {
        $chains[$_] = 80 * $a2[$_] / $grid{num_per_mile}[$_]
          ;    #num_per_mile: even 1 or 1/6 etc. is OK
        if ( my $z = $chains[$_] - int $chains[$_] ) {
            my $g = $a2[$_] % 100;
            my $i = $grid{num_per_mile}[$_] / 80;
            my @j;
            for my $c ( 0 .. 9999 ) {
                push @j, $i * $c;
                last if abs $_[$_] < $j[-1];
            }
            my $k = qq(@j[0..4]) . " .. " . qq(@j[-4.. -1]);
            my $w = <<EOF;

Take a deep breath. This might not be comprehensible, even to me...
You wanted position @_, which we need to use @a2
from where you told us to count from, (@{$grid{origin}{address}}).

$a2[$_] cannot be expressed in an integer of chains (1/80 of a mile)
which is the limit of precision of our program and in fact the PID
_xxxyyy format. Seen with @_; Residue=$z.

Let's see, at $grid{num_per_mile}[$_] numbers per mile,
that means an address points only at
(+-)$k. which won't exactly hit $_[$_]!
If it does then beware of some cross sector problem...
Or something else in this message I haven't refined.
Got it, Holmes? Not me :-(
EOF
## One day show the last few items before the target number too.
            die "$w\n ";    #to still get line number, and also newline.
        }
    }
    my $r = Twsp2mi::point_id_plus_chains $grid{origin}{id}, @chains;
    warn "Failure with @_" unless $r;
    return $r;
}

=head1 EWNS

A suffix maker....

 for my $x ( -1, 0 ) {
     for my $y ( -1, 0 ) {
	 printf "%2d %2d %s\n", $x, $y, EWNS( $x, $y );
     }
 }
 -1 -1 S
 -1  0 N
  0 -1 W
  0  0 E

=cut

sub EWNS {
    substr "EWNS", $_[0] * 2 + ( $_[1] < 0 ), 1;
}

=head1 EWNS2

Takes a pair of numbers, and returns the EWNS versions.

 print join ",", EWNS2(-150,330);
 150,W,330,N

=cut

sub EWNS2 {
    my @ans;
    for my $n ( 0, 1 ) {
        push @ans, abs( $_[$n] ), EWNS( $n, $_[$n] );
    }
    return @ans;
}

=head1 canon

Make address N E S W + numbers canonical, into just (-)number

 600W, 600 W, S 600 all become -600
 600E, 600 E, N 600, 600 all become just 600.

=cut

sub canon {
    my $num;
    for ( $_[0] ) {
        die "$_: no number found" unless /\d+/;
        $num = $&;
        $num *= -1 if /[SW]/;    #No, we're not talking about Franklin St. SW
    }
    return $num;
}

=head1 id2governing_id();

The (single) governing id we define as the id of the (southwest corner of the) section the id belongs to.

Real simple:

 print id2governing_id("AZ230490N0130W0_240433");
 AZ230490N0130W0_200400

=head2 Points lying right upon the north or east edges of sections

It is true that indeed if our target is right on the east or north
edge of a section, we don't really need the neighboring section's
corners. But just to make our programs simpler, we will want them
anyway!

It all has to do with the _xxxyyy structure, and us referring to the
governing section as _x00y00.

=cut

sub id2governing_id {
    my $s = $_[0];
    $s =~ s/_(\d)\d\d(\d)\d\d$/_${1}00${2}00/ or die "Bad args \"@_\"";
    return $s;
}

=head1 surrounding_ids();

The surrounding ids we define as the ids of the four corners of the section the id belongs to.

 print "$_\n" for surrounding_ids("AZ230490N0130W0_240433");
 AZ230490N0130W0_200400
 AZ230490N0130W0_300400
 AZ230490N0130W0_200500
 AZ230490N0130W0_300500

=cut

sub surrounding_ids {
    $_[0] =~ /^(.*_)(\d)\d\d(\d)\d\d$/ or die;
    my @r;
    for my $x ( 0, 1 ) {
        for my $y ( 0, 1 ) {
            push @r, sprintf "%s%d00%d00", $1, $2 + $x, $3 + $y;
        }
    }
    return @r;
}
1;
