#!/bin/bash
#Las Animas Co. CO, but only between PLSS correction lines.
#Needs tons more work.

set -ue
# Actually over such long distances, I should make my lines strings of points with one
#point per intersection... or drop the lines altogether and just make points. Indeed that
#could have symbol: crosshairs.

#Also never mind +proj=pipeline... many things don't support it!

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