#!/usr/bin/perl
# Las Animas Co. CO, but only between PLSS correction lines.
# Here we recognize that Range lines are parallel to meridians of longitude,
# and Township lines are parallel to parallels of latitude.
# Therefore we will no longer do any map projections!
#
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.html
# Created: 2023-11-11T03:22:37+0000
# Last-Updated: 2023-11-14T05:34:52+0000
#     Update #: 98
#

use strict;
use warnings "all";

# OK, let's start from the corner of the
# 7th Guide Meridian West, (R 56/57 W) and
# 6th Standard Parallel South (T 30/31 S).
my %P7GW6PS = ( x => -103.7353009, y => 37.3798999 );

# We must stay between the standard parallels (correction lines).
# Townships: T26-30S. And county extent: R51-67W.

#Populate lookup table
my %adr = ( 0 => { x => ( address("R57W") )[0], y => ( address("T31S") )[0] } );

# It is five townships (30 miles) to the 5th Standard Parallel South.
my %center = %P7GW6PS;
my %range  = ( x => [ -11 * 6 .. 6 * 6 ], y => [ 0 .. 5 * 6 ] );

#use Data::Dumper;print Dumper \%range;exit;
my $increment        = 6;
my $numbers_per_mile = 1000;
my %radius           = ( a => 6378206.4, b => 6356583.8 );    # Clarke 1866

# From analyzing the maps:
sub address {    #uses the "beginning edge" of a township or range.
    my %origin = ( T => -42, R => -70 );
    my $n;
    for (@_) {
        /^(T)(\d+)(S)$/
          || /^(R)(\d+)(W)$/
          || die "\"$_[0]\" should be like T57S, R37W.";
        $n = $2;
        if ( $3 eq "S" || $3 eq "W" ) { $n = -$n; }
        my $adr = -6000 * ( $origin{$1} - $n );
        my $des;
        if    ( $1 eq "T" ) { $des = "N" }
        elsif ( $1 eq "R" ) { $des = "E" }
        else                { die "$1 not T and not R?" }
        return $adr, $des;
    }
}
use Math::Trig;
use Math::Trig ':pi';
sub dlon { 360 * 1609.344 / ( cos deg2rad(pop) ) / $radius{a} / pi2; }
sub dlat { 360 * 1609.344 / $radius{b} / pi2; }

print "WKT,Name\n";
for my $y ( @{ $range{y} } ) {
    for my $x ( @{ $range{x} } ) {
        if ( $y % $increment || $x % $increment ) { next; }
        my %nos;
        my %v = ( x => $x, y => $y );
        for (qw /x y/) {
            $nos{adr}{$_} = $adr{0}{$_} + $numbers_per_mile * $v{$_};
            $nos{cr}{$_}  = ( $_ eq "y" ? -1 : 0 ) + $nos{adr}{$_} / 500;
            die "leftover in CR $_" unless $nos{cr}{$_} = int $nos{cr}{$_};
        }
        printf qq{"POINT (%f %f)",%d/%d\n},
          $center{x} - dlon($y) * $x,
          $center{y} + dlat() * $y,
          $nos{cr}{x}, $nos{cr}{y};
        next;
        printf qq{"POINT (%f %f)",%dE/%d %dN/%d\n},
          $center{x} - dlon($y) * $x,
          $center{y} + dlat() * $y,
          $nos{adr}{x}, $nos{cr}{x}, $nos{adr}{y}, $nos{cr}{y};
    }
}

=x
LANN=37.6442082\ -104.0580960 #Las Animas Co. north notch
#LANN18s=37.3804956\ -104.0631667 #18 Miles south,
# at PLSS corner T30/31S, R59/60W elev 5247 ft., north side of PLSS correction line
#SE https://www.openstreetmap.org/#map=16/37.3818/-103.0748
#NC https://www.openstreetmap.org/#map=16/37.6442/-104.0581
#So SC is https://www.openstreetmap.org/#map=16/37.3818/-104.0581
LABeBa=37.6435097\ -103.0758675 #tri county
center=$LANN pointB=$LABeBa
set $center $pointB
set -- $(geod -p +ellps=GRS80 +units=mi -I -f %.6f <<<$@)
#units: just to see if we got 54 mi = 9 townships x 6 miles...
azimuth=$1 #(perl -wle "print $1 - 180;")
pr_affine=$(
    perl -MMath::Trig -wle '
    $a = deg2rad '$azimuth';
    $c = cos $a;
    $s = sin $a;
    print "+proj=affine +s11=$c +s12=-$s +s21=$s +s22=$c";
    ')
NBSE=37.3818\ -103.0748
#So the numbers are 1000 to the mile, but twice the road names.
#export origin_addresses=60000 27000#Spot chosen by me
origin_addresses=60000\ 90000 #Spot chosen by me
t=/tmp
#nums per mile X miles per m
num_per_meter=$(perl -wle "print 1000 / 1609.344")
CR217_A_14=37.6435061\ -103.1672878
p=placemarks
f=$t/$p
set -- -f CSV -lco GEOMETRY=AS_WKT -oo KEEP_GEOM_COLUMNS=NO
o=$@
set -- $center $origin_addresses
set -- +proj=lcc +lon_0=$2 +lat_0=$1 +lat_1=$1 +x_0=$3 +y_0=$4 \
    +k_0=$num_per_meter
projs=$@ bdl=$rel_id numbers_per_line=1000
bdf=$bdl.kml
#gdalsrsinfo "${projs[0]}" -o wkt2; echo 0 0|proj -v ${projs[0]};exit
echo WKT,Name > $f.csv
for i in 0 1
do
    nohead=1 ../gratWKT 1000 76000 72000 111000 $numbers_per_line $i|
	./road_numbers
done >> $f.csv #bad approach, later jumbled into same KML folder/document
set -- -f LIBKML
kml_opts="$@"
#Use SELECT to reduce clutter and get clean KML!
ogr2ogr $o -s_srs "$projs" -t_srs WGS84 $f.1.csv $f.csv
ogr2ogr $o $f.2.csv $f.1.csv -clipsrc $bdf
for point in LANN LABeBa NBSE CR217_A_14
do
    eval set -- \$$point
    set -- $2 $1
    echo "\"POINT ($@)\",$point"
done >> $f.2.csv #bad approach, later jumbled into same KML folder/document
bdy_as_line=$t/bdy_ln.csv #else opaque white in Google Earth Web...
ogr2ogr $o $bdy_as_line $bdf -nlt LINESTRING
#couldn't run the -sql in ogrmerge.py, so:
s=$t/super.kml
ogr2ogr $kml_opts $s $f.2.csv -sql "SELECT Name FROM \"$p.2\""
set -- $center
viking --map 13 --external $s #--latitude $1 --longitude $2 --zoom 14 $s
#viking --map 13 --external $s # --latitude $1 --longitude $2 --zoom 14 $s
# https://ideditor-release.netlify.app/#background=USGS-Scanned_Topographic&disable_features=boundaries&map=9.74/37.9002/-103.0157