#!/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-02-07T12:17:40+0000
#     Update #: 407

=pod

=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

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 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 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 items.\n",
      scalar keys %sections;    #So one can check if anything's missing...
    return %sections;
}

=head1 parse_seq

For GeoJSONSeq, as made by ogr2ogr.

=cut

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

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 $_;
        use Data::Dump q!pp!;
##	die pp $ph;
        ##        my $name  = $ph->{properties}{FRSTDIVID};
        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] ) }

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...

=cut
