#!/bin/bash
# Highland Park, IL, USA House addressing grids.
# Programs from https://proj.org/ , https://www.generic-mapping-tools.org/
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.html
# Created: 2023-06-25T09:20:23+0000
# Last-Updated: 2023-09-03T10:51:00+0000
#     Update #: 198
set -ue
# https://commons.wikimedia.org/wiki/File:HP_IL_Fire_District_Map_and_oblique_axes_address_grid.jpg
# https://library.municode.com/il/highland_park/codes/code_of_ordinances?nodeId=COOR_TITIXGERE_CH93STSI_ARTVSYSTNU
# Sec. 93.355. - Base lines located. # Sec. 93.360. - System of numbering.
sw=(-87.8210725 42.1525805) #SW corner of city limit
es32=(-87.7700587 42.1523468) #County Line Rd. & Lakeside Place
east_az=$(echo ${sw[@]} ${es32[@]} | gmt mapproject -AF+v -fg -o4)
set -- $(echo ${es32[@]} $east_az 3200 | gmt vector -Ttf -je)
centers=("$*" "$*")
azimuth[1]=$(echo ${sw[@]} $@ | gmt mapproject -AF+v -fg -o4 | gmt convert -i+o-90) #point up
ne=($(echo -87.8083059 42.2034905 $east_az 3000 | gmt vector -Ttf -je)) #From NE 30 point
unset es32 east_az
azimuth[0]=$(echo ${centers[1]} ${ne[@]} | gmt mapproject -AF+v -fg -o4 | gmt convert -i+o0) #already points up
set -- $(gmt math -Q 100 660 DIV 0.3048 DIV =)
NPM=$@
#numbers_per_meter=($@ $@)
numbers_per_meter=( #includes tweaks
    $@ #    $(gmt math -Q -705 -698 DIV $@ MUL =)
    $(gmt math -Q 3600 3614 DIV $@ MUL =) #got the graticule right on the PLSS half section line roads!
)
offset_numbers=(20 0) #seems to make it fit better with the ground truth
#One day, see how different numbers would be if the base lines were in the northwest corner.
for i in 0 1
do
    center=(${centers[$i]})
    set -- +proj=omerc \
	+lonc=${center[0]} +lat_0=${center[1]} \
	+x_0=${offset_numbers[0]} +y_0=${offset_numbers[1]} \
	+alpha=${azimuth[$i]} +gamma=0 \
	+k_0=${numbers_per_meter[$i]}
    projs[$i]=$@
done
if false
then
    cat <<EOF > /tmp/colz
-87.7826490 42.1654800 Roger Williams Avenue & St Johns Avenue
-87.84358 42.194698 1900 W 2300 N
-87.8279551 42.2180667 Old Elm 3600N
-87.78242 42.16613 469 W 750 N
-87.804546 42.183675 800 Central Ave., Moraine Twsp office
-87.794542 42.181368 500 Ravine Dr.
-87.833228 42.180711 1901 York Ln,
-87.817887 42.15627 2000 Old Briar Rd.
-87.815876 42.153782 2000 Bartram Ln.
-87.774535 42.165374 200 Roger Williams Av.
-87.782034 42.171813 300 Cedar Ave.
EOF
    echo E, N: 1>&2
    for i in 0 1
    do
	gmt mapproject -J"${projs[$i]}" -o$i /tmp/colz > /tmp/col$i
    done
    pr -T -m -w 1111 /tmp/col?|perl -anwle \
       'printf "%5.0f %5.0f %s\n", @F[ 0, 1 ], join " ", @F[ 4 .. $#F ];' 1>&2
fi
numbers_per_graticule=(100 100)
#number_range=("-3500 0" "0 3900")
number_range=("-5000 0" "0 5000")
for i in 0 1
do
    ./wkt.pl ${number_range[@]} ${numbers_per_graticule[$i]} $i|
	gmt mapproject -I -J"${projs[$i]}"
done |
    gmt spatial -T$bdyfile |
    gmt 2kml -K -Fl -Wthick,red
gmt 2kml -K -O -Fl -W5p,blue $bdyfile

if : 7/7 8/8 ...
then
    set -x -- ${centers[0]}
    az2=$(gmt math -Q ${azimuth[@]} 270 ADD ADD 2 DIV =)
    N1X=$(gmt math -Q 8 $NPM DIV SQRT =)
    gmt math -T0/2300/${numbers_per_graticule[0]} T $N1X MUL = |
	perl -anwle '@K=('$1','$2','$az2',reverse @F); print "@K";' > /tmp/h
	gmt vector -Tte -je /tmp/h > /tmp/i
	gmt convert /tmp/i /tmp/h -A -o0,1,6|
	gmt 2kml -K -O -N2 -Fs
#Gosh.... I'll never get the above right. And I need to make a universal ADDR=>LON LAT converter.
fi
if :
then
    {
	for i in sw ne
	do
	    eval echo \${$i[@]} $i
	done
	for i in 0 1
	do echo ${centers[$i]} o$i
	done
    } |
	tee /tmp/yy |
	gmt 2kml -O -Nt
fi
#set|grep '^[a-z]'
#gpxsee shows color, viking doesn't, but can be run via shell.