#!/usr/bin/perl
# PLSS point ids to addresses, for Evanston IL
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.htm
# Created: 2024-02-18T01:36:52+0000
# Last-Updated: 2024-03-11T03:28:29+0000
#     Update #: 720
#

=pod

Actually why am I still forcing mile/80 granularity when it comes to giving address points?
Any integer should be fine...

=head1 "Off the deep end" with gdaltransform

Holy moly, now going off the deep end, with gdaltransform. Our
concept today is only give gdaltransform only the four GCPs relevant
to the point in question (the place where we want to put our address
tick.) We give it only the four surrounding section corners.

=cut

use strict;
use warnings q(all);
use PointId2Address;
use Twsp2mi;
use Data::Dumper;
$Data::Dumper::Deepcopy++;
$Data::Dumper::Sortkeys++;
my %co;

=head1 Evanston actually has three or more grids

(if one looks at it from the PLSS.)

=over

=item T41N, south of Central St.

Lots of Evanston. Chopped of before a full mile, at Central St.

=item T42N, which we will define as Central St. and northwards

Note the PLSS data has nicely created missing sections out of the
Ouilmette Reservation, which we will thankfully use.

=item Downtown

Downtown is a part of T41N that we will apply a tilt to...

=back

=cut

my $false_corners;
if ( $ARGV[0] eq '--false-corners' ) {
    $false_corners++;
    shift;
}
my %grids;
$grids{S} = {
    origin => {
        id      => 'IL030410N0140E0_100240',   # corner of Howard & Asbury
        address => [ -1300, 100 ]              # just a random observation point
    },
    num_per_mile => [ 1200, 800 ]              # from City ordinance
};
%{ $grids{N} } = %{ $grids{S} };               #deep copy

# Must start using custom assignments due to distortion, visible on
# the grids on the sides of the official city maps too

# Highland is 3200w, so section line should be 3800w.
# IL030420N0130E0_300100 W city limit 3800 W 2600 N
# IL030420N0130E0_400100 Ewing 2800 W 2600 N
# IL030420N0130E0_500100 Ashland 1500 W 2600 N

$grids{N}{origin} = {
    id      => 'IL030420N0130E0_500100',    # corner of Central & Ashland
    address => [ -1500, 2600 ]              # just a random observation point
};

sub pid2sector {
    return $_[0] =~ /0420N/ ? 'N' : 'S';    #Add Downtown one day?
}

# Here we go. We are reading in the corners!
while (<>) {
    chomp;
    my @F = split /,/;
    $co{ $F[2] } = [ @F[ 0, 1 ] ];
}

=head1 Irregular shaped sections

To top it all off, the study area (Evanston) includes irregular
sections, ragged on the east at the lakeshore, and chopped on the
north at Central St.

Using any irregular section corners would distort our results!

OK, we will build a plank into the lake of false teeth sections.

Maybe Central St. isn't so bad...

=cut

if ($false_corners) {
    my %new_co;

    # It's the world of "false teeth"! We apply our fixes to the
    # lakefront broken edge sections...
    my @needed = qw/
      IL030410N0140E0_100500
      IL030410N0140E0_100600
      IL030410N0140E0_200500
      IL030410N0140E0_200600
      IL030410N0140E0_300200
      IL030410N0140E0_300300
      IL030410N0140E0_300400
      IL030420N0130E0_600100
      IL030420N0130E0_600200
      /;
    for
      ( # Order matters, build in rows, then in each row by columns eastward.
        sort {
            ( Twsp2mi::point_id2miles($a) )[1]
              <=> ( Twsp2mi::point_id2miles($b) )[1]
              || ( Twsp2mi::point_id2miles($a) )[0]
              <=> ( Twsp2mi::point_id2miles($b) )[0]
        } @needed
      ) {
        if ( $co{$_} ) {
            warn "Already know coordinates for $_. Skipping!";
            next;
        }
        warn "Making supplemental $_";
        my @id;
        $id[0] = $_;
        my @ll;
        my @comm = qw/echo/;
        for ( 1, 2 ) {
            $id[$_] = Twsp2mi::point_id_plus_miles( $id[0], -$_, 0 );
            unless ( $ll[$_] = $co{ $id[$_] } ) {
                $id[$_] =~ s/0140E0_100/0130E0_700/;
                unless ( $ll[$_] = $co{ $id[$_] } ) {
                    die "While computing\n$id[0], no known corner at:\n$id[$_]";
                }
            }
        }
        my @G = qw/| geod +ellps=WGS84 +units=us-mi -f %.6f/;
        for ( 2, 1 ) { push @comm, reverse @{ $ll[$_] } }
        push @comm, @G, "-I";
        my $res = qx(@comm) or die;
        my ( $az, $back_az, $dist ) = split /\s+/, $res;
        @comm = ( "echo", ( reverse @{ $ll[1] } ), $az, 1, @G );
        $res  = qx(@comm) or die;
        ( my ( $lat, $lon ), $back_az ) = split /\s+/, $res;
        @{ $new_co{$_} } = @{ $co{$_} } = ( $lon, $lat );
    }
    for ( keys %new_co ) {
        print join ",", @{ $new_co{$_} }, $_;
        print $/;
    }
    exit;
}
my @targets = generate_sample_ticks();

sub generate_sample_ticks {
    my @ticks;    # Example ticks
    for ( my $i = 400 ; $i <= 2400 ; $i += 300 ) {
        push @ticks, [ $i, -$i, 100 ];    #Howard St.
        ##[LABEL, X, Y]
    }
    for ( my $i = 700 ; $i <= 3600 ; $i += 300 ) {
        push @ticks, [ $i, -$i, 2100 ];    #Simpson St.
    }
    for ( my $i = 600 ; $i <= 3600 ; $i += 300 ) {
        push @ticks, [ $i, -$i, 2600 ];    #Central St.
    }
    for ( my $i = 300 ; $i <= 2599 ; $i += 200 ) {
        push @ticks, [ $i, -1300, $i ];    #Asbury Ave. S sector
    }
    for ( my $i = 2600 ; $i <= 2800 ; $i += 200 ) {
        push @ticks, [ $i, -1305, $i ];    #Asbury Ave. N sector
    }
    for ( my $i = 300 ; $i <= 1100 ; $i += 150 ) {
        my $k = 100 + $i;
        push @ticks, [ $k, -$k, $k ];      #West=North Diagonal
    }
    return @ticks;
}

my %tarr;
for (@targets) {
    my ( $label, @addr ) = @$_;
    my $id =
      PointId2Address::addr2id( $grids{ addr_pair2sector(@addr) }, @addr );
    $tarr{$id}{gov} = PointId2Address::id2governing_id($id);
    unless ( @{ $tarr{$id}{coordinates} } =
        govid_and_adr2coords( $tarr{$id}{gov}, @addr ) ) {
        warn "Skipping $label @addr";
        next;
    }
    die "Cannot make $id coordinates" unless 2 == @{ $tarr{$id}{coordinates} };
    print join ",", @{ $tarr{$id}{coordinates} }, $label;
    print $/;
}

sub govid_and_adr2coords {
    my ( $govid, @adr ) = @_;
    my @surrounding_ids = PointId2Address::surrounding_ids($govid);
    my @args;
    for (@surrounding_ids) {
        push @args, "-gcp";
        push @args, id2addr_local($_);
        unless ( $co{$_} ) {
            my $mess = "@args";
            for ($mess) {
                s/-gcp/\n$&/g;
                s/\s+\n/\n/g;
                s/^\n//;
            }
            warn "Gov $govid, adr @adr:\n$mess: No known corner of $_";
            return undef;
        }
        push @args, @{ $co{$_} };
    }
    my $result = qx!echo @adr | gdaltransform -output_xy @args!;
    die unless $result;    #poor error checking
    chomp $result;
    return split / /, $result;
}

sub id2addr_local {
    my $id     = $_[0];
    my $sector = pid2sector($id);
    ##Hmm, accessing a global variable...
    $grids{$sector}{target} =
      { id => $id };       # Must wipe out previous {target}{miles}
    my @o = PointId2Address::id2addr_raw( $grids{$sector} ) or die;
    if ( $id =~ /600$/ ) {    #Central street!
        $o[1] = 2600;
    }
    return @o;
}

=head1 addr_pair2sector() - determine via address what sector of Evanston this is.

Lets call our sectors S, for south of Central St., N for Central St. and northward, and D for Downtown.

=cut

sub addr_pair2sector {
    my ( $x, $y, @rest ) = @_;
    if ( $y >= 2600 ) {
        return 'N';
    }
    elsif (0) { return 'D' }    #not ready yet
    else      { return 'S' }
}