#!/usr/bin/perl
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.html
# Created: 2023-12-12T22:14:09+0000
# Last-Updated: 2024-05-21T01:55:04+0000
#     Update #: 586

=head1 NAME

PLSS_corners_4 - Find just four corners of US PLSS sections, without needing centers

=head1 SYNOPSIS

See PLSS_corners.

=head2 Input

Automatically detected, from among,

=over

=item *

BLM JSON... the same format as used in PLSS_corners.pm .

=item *

GeoJSON e.g., as made by ogr2ogr.

=back

=head1 DESCRIPTION

Given a ring of nodes representing a US PLSS section, find the four corner nodes.

=head1 METHOD 1

We first make a bounding box, then compute which nodes are closest to the box corners.
These must be the four corners of the PLSS section.

=head1 METHOD 2

(sniff4()): Better yet, instead We make a table of each node, and the rings it belongs to.

The nodes that belong to four rings are the ones we are pursuing.

Well, the ID of the northwest ring's southwest corner is the ID for this point!

=head1 Filters

See the code for various optional environment variables... which I am
removing because I can't tell what I meant.

=cut

package PLSS_corners_4;
use strict;
use warnings q(all);
use JSON::PP;             # Installed already...
use Math::Trig;
use Math::Trig qw(:pi :great_circle deg2rad rad2deg);
use PLSS_corners;
use Twsp2mi;
use Data::Dump qw.pp.;    #Uh oh, non necessarily installed already...
sub NESW { deg2rad( $_[0] ), deg2rad( 90 - $_[1] ) }

sub parse {
    my $ph;
    { local $/; $ph = decode_json <>; }
    my $count;
    my %sections;
    ##We sniff the first node to probe contents format,
    if ( $ph->{features}[0]{geometry}{rings} ) {
        ## Parse BLM PLSS JSON
        for ( @{ $ph->{features} } ) {
            my $name = $_->{attributes}{landdescription};
            $sections{$name} = $_->{geometry}{rings}[0];
            $count++;
            if ( $ENV{limit} ) { last if $count >= $ENV{limit}; }
        }
    }
    elsif ( $ph->{features}[0]{geometry}{coordinates} ) {
        ## parse GeoJSON, e.g., that was made by ogr2ogr
        for ( @{ $ph->{features} } ) {
            my $name = $_->{properties}{FRSTDIVID};
            $sections{$name} = $_->{geometry}{coordinates}[0][0];
            if ( $ENV{limit} ) { last if $count >= $ENV{limit}; }
        }
    }
    else {
        die
          "Seems like that flavor of JSON hasn't been added to my program yet.";
    }
    printf STDERR "Found %d sections.\n",
      scalar keys %sections;    #So one can check if anything's missing...
    return %sections;
}

my %sides;
for ( $ENV{sides} ) {
    next unless length;
    for ( split /,/ ) {
        $sides{$_}++;
    }
}

=head1 parse_seq

For GeoJSONSeq, as made by ogr2ogr.

=cut

sub parse_seq {
    if ( $ENV{intervals} ) {
        my $m = 80 / $ENV{intervals};
        if ( $m != int $m ) {
            my $mes = <<EOF;
We've extended the PLSS Points ids, so that instead of just expressing
down to 1/16 of a mile, _x05x05, we go down to 1/80, _x01x01 .

However you Sir have used intervals=$ENV{intervals} and 80 /
intervals=$ENV{intervals} is not an integer. So we can't deal cleanly
with that.
EOF
            die substr $mes, 0, -1;    #Hack to force line number.
        }
        die "When using intervals= also use sides=1 sides=1,3 etc."
          unless %sides;
    }
    while (<>) {
        my $ph   = decode_json $_;
        my $name = $ph->{properties}{FRSTDIVID}
          || $ph->{properties}{PLSSID}    #Township mode
          ;
        my @nodes = @{ $ph->{geometry}{coordinates}[0] };
        pop @nodes;                       #remove repeating closing node
        for my $c ( corners( $name, @nodes, ) ) {
            printf "%s,%s,%s\n", @{ $c->[1] }, $c->[0];
        }
    }
}

sub parse_seq_MN46 {
    ## Hmm, we encountered a different Minnesota format from
    ## https://gisdata.mn.gov/dataset/plan-mngeo-trs . I guess I need to
    ## make a special function... hmm, the nation probably has many
    ## special formats...
    while (<>) {
        my $ph = decode_json $_;
        use Data::Dump q!pp!;
        my $name = sprintf "%s%02d%03d0N%03d0W0SN%02d0", "MN", 46,
          @{ $ph->{properties} }{qw/TOWN RANG SECT/};
        my @nodes = @{ $ph->{geometry}{coordinates}[0] };
        pop @nodes;    #remove repeating closing node
        for my $c ( corners( $name, @nodes, ) ) {
            printf "%s,%s,%s\n", @{ $c->[1] }, $c->[0];
        }
    }
}

sub go {

    # warn "This part of the program was just an early stage...";
    # Actually only tested on townships...
    my %sections = parse();
    for ( keys %sections ) {
        my @nodes = @{ $sections{$_} };
        pop @nodes;    #remove repeating closing node
        for my $c ( corners( $_, @nodes, ) ) {
            printf "%s,%s,%s\n", @{ $c->[1] }, $c->[0];
        }
    }
}

sub corners {
    my ( $name, @nodes, ) = @_;
    my @zz;
    for ( 0, 1, ) {
        my @ranked = sort { $a->[$_] <=> $b->[$_] } @nodes;
        push @zz, $ranked[0][$_], $ranked[-1][$_];
    }
    my ( $w, $e, $s, $n, ) = @zz;    #trying to keep the code readable.
    my @extremes = ( [ $e, $n, ], [ $e, $s, ], [ $w, $s, ], [ $w, $n, ], );
    my @corners;
    for ( 0 .. 3 ) {
        my @t       = NESW( @{ $extremes[$_] } );
        my @nearest = sort {
            great_circle_distance( @t, NESW(@$a), )
              <=> great_circle_distance( @t, NESW(@$b), )
        } @nodes;
        ##could have remove the node after we confirm it was a corner.
        ##could do NESW much earlier to avoid repeating-- if it matters.
        $corners[$_] = [ PLSS_corners::point_id( $name, $_, ), $nearest[0], ];
    }
    if ( $ENV{intervals} ) {
        for my $q ( 0 .. $#corners ) {
            if ( $sides{$q} ) {
                my @waypoints =
                  interpolate( $ENV{intervals},
                    $corners[$q][1], $corners[ ( $q + 1 ) % 4 ][1],
                  );
                for (@waypoints) {
                    push @corners,
                      [
                        PLSS_corners::point_id(
                            $name, $q, -80 * $_->[0] / $ENV{intervals},
                        ),
                        $_->[1],
                      ];
                }
            }
        }
    }
    return @corners;
}

=head2 interpolate

Given the number of slices we want, and two end points
return each slice # and intervening point.

=cut

##Don't complicate this with any more parameters!
sub interpolate {
    use Math::Trig 'great_circle_waypoint';
    my ( $intervals, @ends_deg, ) = @_;
    my @ends;
    for (@ends_deg) { push @ends, NESW @$_; }
    my @waypoints;
    for ( 1 .. $intervals - 1 ) {
        push @waypoints,
          [ $_, [ SWNE( great_circle_waypoint( @ends, $_ / $intervals, ) ) ], ];
    }
    return @waypoints;
}

sub rad2deg2 { sprintf "%.7f", rad2deg $_[0]; }
sub SWNE { rad2deg2( $_[0] ), 90 - rad2deg2( $_[1] ) }

=head1 Township mode

Here we feed it township rings instead of section rings.

On output we assume the four extreme corners.

 ogr2ogr -f GeoJSONSeq big.json -sql \
 "SELECT PLSSID FROM PLSSTownship LIMIT 3" IL_CadNSDI_V2.gdb.zip
 perl -MPLSS_corners_4 -we 'PLSS_corners_4::parse_seq;' big.json
 -89.972452,41.4962174,IL040170N0040E0_700700
 -89.9723187,41.4095994,IL040170N0040E0_700100
 -90.0864564,41.4106578,IL040170N0040E0_100100%
 -90.0849627,41.4977094,IL040170N0040E0_100700...

=head1 sniff4()

This is where we don't use any Math::Trig at all, but instead use an "AI" approach....
mentioned at the top of this page.

=cut

sub sniff4 {
    die if @_;
    my %sections = parse();
    my %points;
    my @confirmed;
    my %tally;
    for my $sid ( keys %sections ) {
        my @ring = @{ $sections{$sid} };
        pop @ring;    #drop closing node
        for (@ring) {
            push @{ $points{qq(@$_)} }, $sid;
        }
    }
    for ( sort keys %points ) {    #sort: for consistent bug label numbers
        my $touches = @{ $points{$_} };
        $tally{$touches}{total}++;
        if ( $touches == 4 or $touches == 3 ) {    #Main case
            next if $ENV{debug} && $touches == 4;
            next if $touches == 4;
            my $pid = northeast( @{ $points{$_} }, $tally{$touches}{total}, );
            ## When debugging make the good ones have no label for easier viewing of the junk ones
            if ($pid) {
                if ( $ENV{debug} ) { $pid = "" }
            }
            else {
                next unless $ENV{debug};
                $pid = $tally{$touches}{total};
            }
            ## Junk caused by four corner points at:

            ## Section corners along E-W correction lines that just happen to
            ## be at exact jog lengths of one mile, two miles, etc.

            ## Section corners along N-S meridian zone boundaries
            ## (e.g., IL02/IL03) where the four sections just
            ## happen to meet.

            ## Other cases.

            ## Sure, we could tweak our program one day to salvage
            ## these often close to 640 acre sections one day...
            ## if it was worth the trouble. We could indeed even
            ## check GISACRES etc. fields, if we really want to
            ## screw things up... or OGR ST_Area(GEOMETRY)...
            ## or https://gdal.org/user/ogr_sql_dialect.html#ogr-geom-area

            if ( $pid ne 0 ) {
                my @l = split / /, $_;
                push @confirmed, [ @l, $pid, ];
                $tally{$touches}{accepted}++;
            }
        }
        elsif ( $touches == 3 )
        {    #Still very valuable points of correction lines, etc.
            my $pid = best_of_3( @{ $points{$_} }, $tally{$touches}{total}, );
            if ( $pid ne 0 ) {
                my @l = split / /, $_;
                push @confirmed, [ @l, $pid, ];
                $tally{$touches}{accepted}++;
            }
        }
        elsif ( $touches == 2 )
        { #Run of the mill points on section lines, but not at a corner. Discard.
        }
        elsif ( $touches == 1 )
        {    #Rare, found only on the outer corners of our sample area. Discard.
        }
        else {
            die "$_: same point on edge of $touches different sections?";
        }
    }
    print STDERR "Corners accepted/total\n";
    for ( reverse sort keys %tally ) {
        printf STDERR "$_: %d/%d\n", $tally{$_}{accepted} || 0,
          $tally{$_}{total} || 0;
    }
    use Data::Dumper;
    $Data::Dumper::Sortkeys++;
    return @confirmed;
}

=head1 sniff4p()

Like sniff4p but prints CSV rows.

=cut

sub sniff4p {
    my @o = sniff4(@_);
    local $, = ",";
    for (@o) {
        print @$_;
        print $/;
    }
}

=head1 northeast()

Given three or four adjacent section SNids and a debugging serial
number (for putting in the map labels when debugging), return the PID
of the southwest corner of the northeast one.

First we make sure they are all using the same meridian.

Then we compute how far they are from that meridian's origin. Then we
inch along them until we are sure we are on the northeast one...

We return the pid of the most northest section, or 0 for failure.

=cut

sub northeast {
    ## Working great so far for 4 corners.
    ## Needs to be made ready for three.

    use Data::Dumper;
    my $serial_number = pop;
    unless ( same_meridian(@_) ) {
        warn "Marker $serial_number: different meridians";

        # We are only interested in regularly spaced grid points a
        # mile apart. One day in fact could salvage some of these.
        return 0;
    }
    for (@_) {
        if ( non_regular_section_id($_) ) {
            warn "Marker $serial_number: non_regular_section_id $_";
            return 0;
        }
    }
    my %distance;
    for (@_) {
        my $pid = Bou2colrow::trsec_dir2id( $_, q(SW), );
        @{ $distance{$pid} } = Twsp2mi::point_id2chains($pid);
    }
    my @sorted =
      sort {    # First sort n/s, then e/w
             $distance{$a}->[1] <=> $distance{$b}->[1]
          || $distance{$a}->[0] <=> $distance{$b}->[0]
      } keys %distance;
    ## Now we start checking. We expect SW section to NE sections have 80, 80 chain differences.
    my @perfect;
    for ( 0, 1, ) {
        $perfect[$_] =
          0 +
          ( ( $distance{ $sorted[-1] }[$_] - $distance{ $sorted[0] }[$_] ) ==
              80 );
    }
    if ( $perfect[0]  && $perfect[1] ) { return $sorted[-1]; }
    if ( !$perfect[0] && !$perfect[1] ) {
        my $m = <<~EOF;
        Marker $serial_number: ids not 1 mi E nor N.
        Half the time the the BLM recorded the numbers off the topo
        map wrong!
        EOF
        $m .= present(@_);
        chomp $m;
        warn $m;
        return 0;
    }
    if ( neighboring_n_s_townships(@sorted) ) { return $sorted[-1] }
    my $m = <<~EOF;
    Marker $serial_number: Perhaps in fact a BLM error. Check that the
    section numbers they use here are correct.
    EOF
    $m .= present(@_);
    chomp $m;
    warn $m;
    return 0;
}

=head1 same_meridian()

Returns 1 if all the IDs given to it belong to the same meridian.
E.g., [IL02 and IL03] are different. [IL02 and IL02] are the same, and
even [IL02 and IN02] are the same also!

=cut

sub same_meridian {
    ## Yes, we will ignore same or different states! We only examine meridians.
    my %meridian;
    for (@_) { $meridian{ substr( $_, 2, 2, ) }++ }
    return keys %meridian == 1;
}

=head1 non_regular_section_id()

Catch e.g., T032AN, SN35A, etc.

returns 1 if non-regular, else 0.

=cut

sub non_regular_section_id {    # one day accept A etc. instead of 0's.

    return $_[0] !~ /^[A-Z][A-Z]\d{2}\d{3}0[NS]\d{3}0[EW]0SN[0-3]\d0$/;
}

=head1 neighboring_n_s_townships()
is this a case of e.g., (T11N) T11N T12N T12N?

Given three or four ...SN..'s,
If all are in the same TxxN or S: return 0;
If two are north and one or two south, return 1;
If one is north and two are south, return 0;
If any is farther than one twsp away, return 0 and warn.

=cut

sub neighboring_n_s_townships {
    my %chains_north;
    for (@_) {
        $chains_north{ ( Twsp2mi::point_id2chains($_) )[1] }++;
    }
    unless ( keys %chains_north == 2 ) {
        warn "Not two different of townships north, only ",
          scalar keys %chains_north, $/, present(@_), " ";   #to get line number
        return 0;
    }
    ##Only considering six mile granularity, that way if a township is
    ##cut off below the first row of sections, it will still work...
    my @k        = sort keys %chains_north;
    my $distance = $k[1] - $k[0];
    if (   $chains_north{ $k[0] } == 2
        && $chains_north{ $k[1] } == 1 ) {

        #reject, a bottom heavy 700 case we are not interested in this year.
        return 0;
    }
    if ( $distance == 80 * 6 ) {
        return 1;
    }
    else { return 0 }
}

## sub neighboring_n_s_townships {
##     my %townshipsNS;
##     for (@_) {
##         my @extract = Twsp2mi::pid2extract( $_, 1, );    #1: N/S
##         $townshipsNS{qq(@extract[ 0, 1, ] )}{count}++;
##     }
##     unless ( keys %townshipsNS == 2 ) {
##         warn "Not two different of townshipsNS, only ",
##           scalar keys %townshipsNS, $/, present(@_), " ";    #to get line number
##         return 0;
##     }
##     my @u;
##     for ( keys %townshipsNS ) {
##         push @u, [ substr( $_, 0, 3, ), substr( $_, 3, 1, ), 1, 0, 0, ];
##     }
##     ##Only considering six mile granularity, that way if a township is
##     ##cut off below the first row of sections, it will still work...
##     die Dumper \@_, \@u, \%townshipsNS;
##     my $dist      = Twsp2mi::twdiffi_chains(@u);
##     my $neighbors = ( abs $dist == 80 * 6 );
##     unless ($neighbors) {
##         warn 77, $dist, Dumper \@_;
##         return 0;
##     }
##     return $neighbors;
## }

sub present {    #make things present-able.
    my $str;
    $str .= sprintf "$_,$/" for sort @_;
    return $str;
}
1;