#!/usr/bin/perl
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.html
# Created: 2023-11-17T13:45:38+0000
# Last-Updated: 2024-05-16T02:45:54+0000
#     Update #: 904

=head1 NAME

PLSS_corners - Find corners of US PLSS sections

=head1 SYNOPSIS

With standard input being BLM output as seen in
https://www.jidanni.org/geo/house_numbering/grids/utilities/fetch_blm_sections ,

 use PLSS_corners;
 PLSS_corners::process();

Output is CSV: Longitude,Latitude,Name
but without the header.

=head1 DESCRIPTION

Given a ring of nodes representing a US PLSS section, and its center,
compute which nodes are the corners, and which are the half mile
points along the edges.

Compute the corners, of PLSS sections obtained from JSON of
https://gis.blm.gov/arcgis/rest/services/Cadastral/BLM_Natl_PLSS_CadNSDI/MapServer/exts/CadastralSpecialServices/GetLatLon

Because that site has lots of section data, but we couldn't find a
nationwide corner retrieving website like it.

=head1 METHOD

For each quadrant, the longest distance from the center of the
section (given in the data) to a node must be to a (reasonable choice
for) corner node for that quadrant.

Indeed, we can do the same method, using quadrants in an "x" pattern,
instead of "+", to get the half mile nodes. These one expects have the
shortest distance from the center. Note not all input data contains
half mile nodes, so the quality of these might vary...

And if we combine "+" and "x", we get a skew compromise that allows us
to use just one bucketing loop, bringing the complexity of our program
back to bearable levels. Four buckets instead of 4 * 2 buckets.

Corner node names: we choose "Bureau of Land Management Geographic
Coordinate Database (GCDB) Publication Point IDs" format, e.g.,
`600100`, referring to the southwest corner of section 36. The half
mile node east of it would be called 640100.

Setting environment variable perfect_only skips all problems sections.

=head1 Quadrants

We refer to them as:

=head2 Corners

=over

=item 0 NE

=item 1 SE

=item 2 SW

=item 3 NW

=back

=head2 Sides

=over

=item 0 E

=item 1 S

=item 2 W

=item 3 N

=back

=cut

package PLSS_corners;
use strict;
use warnings q(all);
use Data::Dumper;
##$Data::Dumper::Indent   = 3;
$Data::Dumper::Sortkeys = 1;
use JSON::PP;    # Installed already...
use Math::Trig;
use Math::Trig qw(:pi :great_circle deg2rad rad2deg);
use Bou2colrow;
sub NESW { deg2rad( $_[0] ), deg2rad( 90 - $_[1] ) }
my $skew = atan( 1 / 2 );

=head1 Bucket details

We place our bucket boundaries at spots 200160 and every 90 degree
rotation: 160100, 100120, 120200, using Sec 31 as an example.
Brilliantly including just one corner and one half mile point in each
bucket! (Bucket boundaries at 26 deg. azimuth (where we cross from
bucket #3 to #0) and every 90 degree increment.

=cut

sub buck { ( ( $_[0] - $skew + pi2 ) / pip2 ) % 4 }
my $slack = deg2rad 3;

sub check_angles {
    my @t;
    for (@_) {    #solved LAST NODE loop back
        next unless $_->{label};
        push @t, $_;
    }
    if ( @t % 4 ) {
        return 1;
    }
    my $previous = $t[-1];
    for (@t) {
        return 1
          if $_->{azimuth}{radians} -
          $previous->{azimuth}{radians} -
          ( @t == 8 ? pip4 : pip2 ) >= $slack;
        $previous = $_;
    }
    return 0;
}

sub process {
    ## Doesn't work for townships yet, only sections. So please don't
    ## feed it townships, else results are a mess. See PLSS_corners_4
    ## instead.
    my $ph;
    my @sections;
    { local $/; $ph = decode_json <>; }

    for ( 0 .. $#{ $ph->{coordinates} } ) {
        if ( $ENV{avoid} ) {    #a pattern of plssids to skip
            next if $ph->{coordinates}[$_]{plssid} =~ $ENV{avoid};
        }
        if ( $ENV{filter} ) {    #a single plssid to debug
            next unless $ENV{filter} eq $ph->{coordinates}[$_]{plssid};
        }
        my %cur_sec;
        for my $i (qw/plssid/) {
            $cur_sec{$i} = $ph->{coordinates}[$_]{$i};
        }
        for my $i (qw/lon lat/) {
            push @{ $cur_sec{coordinates}{degrees} },
              $ph->{coordinates}[$_]{$i};
        }
        @{ $cur_sec{coordinates}{radians} } =
          NESW @{ $cur_sec{coordinates}{degrees} };
        for my $c ( 0 .. $#{ $ph->{features}[$_]->{geometry}->{rings}[0] } - 1 )
        { ## -1: last one just closes the loop and is a repeat. I could check though someday.
            push @{ $cur_sec{nodes}[$c]{coordinates}{degrees} },
              @{ @{ $ph->{features}[$_]->{geometry}->{rings}[0] }[$c] };
            @{ $cur_sec{nodes}[$c]{coordinates}{radians} } =
              NESW @{ $cur_sec{nodes}[$c]{coordinates}{degrees} };
        }
        push @sections, \%cur_sec;
    }
  S: for my $s (@sections) {
        my @buckets;
        for my $n ( @{ $s->{nodes} } ) {
            $n->{azimuth}{radians} =    #from center
              great_circle_direction(
                @{ $s->{coordinates}{radians} },
                @{ $n->{coordinates}{radians} }
              );
            $n->{distance} =            #from center, radians
              great_circle_distance(
                @{ $s->{coordinates}{radians} },
                @{ $n->{coordinates}{radians} }
              );
            push @{ $buckets[ buck( $n->{azimuth}{radians} ) ] }, $n;
        }
        for ( 0 .. $#buckets ) {
            my @m;
            if ( $buckets[$_] ) { push @m, 0; }
            else {
                if ( $ENV{perfect_only} ) { next S; }
                ## Garbage in, garbage out.
                mywarn( sprintf "%s: No MILE at %s",
                    $s->{plssid}, point_id( $s->{plssid}, $_, 0 ) );
                next;
            }
            @{ $buckets[$_] } =    #Longest is corner, shortest is side
              sort { $b->{distance} <=> $a->{distance} } @{ $buckets[$_] };

            if ( $buckets[$_][1] ) { push @m, 1; }
            else {
                if ( $ENV{perfect_only} ) { next S; }
                mywarn( sprintf "%s: No half at %s",
                    $s->{plssid}, point_id( $s->{plssid}, $_, 1 ) );
            }
            for my $ls (@m) {
                @{ $buckets[$_] }[ -$ls ]->{label} =
                  point_id( $s->{plssid}, $_, $ls );
            }
        }
        if ( check_angles( @{ $s->{nodes} } ) ) {
            if ( $ENV{perfect_only} ) { next S; }
            mywarn( sprintf "%s: has distortions", $s->{plssid} );
        }

        if ( $ENV{filter} ) {
            if ( $ENV{stop} ) {
                for (@buckets) {    #easier to read:
                    for (@$_) {
                        $_->{azimuth} = rad2deg $_->{azimuth}{radians};
                        delete $_->{coordinates};
                    }
                }
                die Dumper \@buckets;
            }
        }
    }
    for (@sections) {
        for ( @{ $_->{nodes} } ) {
            next unless $_->{label};
            print join ",", @{ $_->{coordinates}{degrees} }, $_->{label};
            print "\n";
        }
    }
}

=head1 CAVEATS

Yes, expect that neighboring sections will produce
2nd copies of shared points. So, the user can run "sort(1) -u"
on them. Or next time just get e.g., sections 33, 35,
if not interested in the half mile points of 34.

Also 700xxx will be a duplicate of 100xxx in the
next township. So use grep -v 700...$ . Yes,
all expected.

=head1 point_id

Given e.g., PLSS id WA330200N0280E0SN340 and quadrant 3 (NW), return
WA330200N0280E0_400200. Or return the following half mile, headed for
the next quadrant corner, if given extra parameter of 1.

 $ perl -MPLSS_corners -lwe '
 for my $q ( 0 .. 3 ) {
     for my $h ( 0, 1 ) {
         print "Q$q H$h ",
           PLSS_corners::point_id( "WA330200N0280E0SN340", $q, $h );
     }
 }'
 Q0 H0 WA330200N0280E0_500200
 Q0 H1 WA330200N0280E0_500140
 Q1 H0 WA330200N0280E0_500100
 Q1 H1 WA330200N0280E0_440100
 Q2 H0 WA330200N0280E0_400100
 Q2 H1 WA330200N0280E0_400140
 Q3 H0 WA330200N0280E0_400200
 Q3 H1 WA330200N0280E0_440200

If that extra parameter is instead negative, it is taken as the 80th of a mile,
yes, also starting at the given quadrant corner and heading to the next, along the edge.

In fact the 1 is just an alias for -40, and all this is for
compatibility with my previous versions...

 perl -MPLSS_corners -we \
 'for my $q(0..3){for my $h(reverse -79,-78,-8){
  printf "Q%d H%3d %s\n", $q, $h,
 PLSS_corners::point_id("WA330200N0280E0SN340",$q, $h);}}'
 Q0 H -8 WA330200N0280E0_500172
 Q0 H-78 WA330200N0280E0_500102
 Q0 H-79 WA330200N0280E0_500101
 Q1 H -8 WA330200N0280E0_472100
 Q1 H-78 WA330200N0280E0_402100
 Q1 H-79 WA330200N0280E0_401100...

Also works for townships,

 for my $q ( 0 .. 3 ) {
     print "Q$q  ", PLSS_corners::point_id( "WA330200N0280E0", $q );
 }

 Q0  WA330200N0280E0_700700
 Q1  WA330200N0280E0_700100
 Q2  WA330200N0280E0_100100
 Q3  WA330200N0280E0_100700

=cut

sub point_id {
    if ( $_[0] =~ /(.*)SN(\d\d).$/ ) {
        my $twsp = $1;
        my $sec  = $2;
        my ( $x, $y ) = Bou2colrow::b2crUS($sec);
        ## Maybe the following could be done much more elegantly with some sine function...
        if ( $_[2] ) {    #fractions of miles
            die "$_[2]: No fractions themselves please..."
              unless $_[2] == int $_[2];
            my $eightieths;
            if    ( $_[2] == 1 )               { $eightieths = 40; }
            elsif ( $_[2] > -80 && $_[2] < 0 ) { $eightieths = abs $_[2]; }
            else { die "$_[2]: Currently only accepting args -79 .. 1"; }
            my @tf;            #Tenths of furlongs! 80 = one mile.
            for ( $_[1] ) {    #quadrant
                if ( $_ == 0 ) {
                    $x++;
                    $tf[1] += ( 80 - $eightieths );
                }              #E quadrant
                elsif ( $_ == 1 ) {
                    $tf[0] += ( 80 - $eightieths );
                }              #S quadrant
                elsif ( $_ == 2 ) { $tf[1] += $eightieths; }         #W quadrant
                elsif ( $_ == 3 ) { $y++; $tf[0] += $eightieths; }   #N quadrant
                else              { die "Quadrant $_?! Use 0..3." }
            }
            return sprintf "%s_%d%02d%d%02d", $twsp, $x,
              ( $tf[0] || 0 ), $y, ( $tf[1] || 0 );
        }
        else {    #maybe this is just a special case of the above...
            for ( $_[1] ) {
                if    ( $_ == 0 ) { $x++; $y++; }    #NE quadrant
                elsif ( $_ == 1 ) { $x++; }          #SE quadrant
                elsif ( $_ == 3 ) { $y++; }          #NW quadrant
            }
            return sprintf "%s_%d00%d00", $twsp, $x, $y;
        }
    }
    elsif ( length $_[0] == 15 ) {    #Township mode
        my ( $x, $y );
        for ( $_[1] ) {
            if    ( $_ == 0 ) { $x = $y    = 7; }    #NE quadrant
            elsif ( $_ == 1 ) { $x = 7; $y = 1; }    #SE quadrant
            elsif ( $_ == 2 ) { $x = $y    = 1; }    #SW quadrant
            elsif ( $_ == 3 ) { $x = 1; $y = 7; }    #NW quadrant
            else              { die "$_[1] supposed to be within 0..3." }
        }
        return sprintf "%s_%d00%d00", $_[0], $x, $y;
    }
    else { die "$_[0]: odd format."; }
}

sub mywarn {
    print STDERR "@_\n"
      ;    #maybe I wanted to add the newline to get rid of the line number...
}

=head1 centers

Given the same json input as process() above,
just return a hash of the PLSS section centers,

 {
    'IL030420N0120E0SN010' => {
                                'lat' => '42.1451893748872',
                                'lon' => '-87.7798724654572'
                              },
    'IL030420N0120E0SN020' => {
                                'lat' => '42.1453073321849',
                                'lon' => '-87.7994410215264'
                              },...

=cut

sub centers {
    ## Only tested on sections, not townships yet.
    my $ph;
    { local $/; $ph = decode_json <>; }
    my %h;
    for ( @{ $ph->{coordinates} } ) {
        for my $l (qw/lat lon/) {
            $h{ $_->{plssid} }{$l} = $_->{$l};
        }
    }
    return %h;
}

1;
