#!/bin/bash
# winnetka.sh --- This file contains some early attempts I used.
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.html
# Created: 2023-09-03T10:54:37+0000
# Last-Updated: 2024-02-17T22:53:53+0000
#     Update #: 7
#
# Even some of the "observations" are less correct that the other
# files in this directory.

set -x
set -ue
s=winnetka

# https://learnosm.org/en/osm-data/getting-data/
# https://codelibrary.amlegal.com/codes/winnetka/latest/winnetka_il/0-0-0-24508
# Village Code
# Section 15.44.010   Numbering buildings.
set -x
#We look at old maps; W border is 1800, O point is 200, 3 miles away
#Right on PLSS section lines.
numbers_per_mile=$(gmt math -Q 1600 3 DIV =)
meters_per_mile=1609.344
set -- $(gmt math -Q $numbers_per_mile $meters_per_mile DIV =)
#in fact it is the same for N/S
numbers_per_meter=($@ $@)

#X observation: on OSM: TWSP 8m VIL 2m 1800W, between 600 to 1000 N.
#Glencoe 36m 1000W, between 1100 and 1400 N.

#y observation: 1000N 30m {Ashbury Ct. and Pine Tree Ln. cul du sacs. middle line}
#Make  500W aligned with Myrtle Av. west line
#Make  700W aligned with Forest Av.
#Make 1000W aligned with Locust Av.
#Make 1400W aligned with Vernon Av.
tweak=(.99 .975)
for i in 0 1
do numbers_per_meter[$i]=$(gmt math -Q ${numbers_per_meter[$i]} ${tweak[$i]} MUL =)
done
offset_numbers=(-200 100)

#2mi e. of Hibbard, 1 mi e. of ridge, origin at turn of winnetka av.

#Origin point via https://maps.cookcountyil.gov/cookviewer
O=(-87.712559 42.093458) #Section corner at big turn ~200 Winnetka Av.
sw=(-87.751560 42.093955) #Hill & Hibbard

place=Winnetka\ IL\ USA

# IT WOULD HAVE BEEN WONDERFUL,
# to use proj's +units.
# But it is not clear,
# how this affects +x_0, +y_0, and +k_0,
# Therefore, I shall refrain from the convience of +units,
# And only feed proj meters, no matter on the command line,
# or on STDIN.

# Also it makes the code much easier to understand! So good. Only use meters.

#We do X and Y independently.
#They might need different +k_0, and they might not be at 90 angles in some cities...
set -eux
# Get the azimuth of our X axis
set -- $(echo ${sw[@]} ${O[@]} | gmt mapproject -AF+v -fg -o4 | gmt convert -i+o-90)
azimuth=($@ $@) #both the same for this town

#We don't care about the 1000 ft, 2000 ft...
#We just observe our origin point is at (100 -200)

#I bet the 1799 etc. Green Bay Road are from the Winnetka system.
D=/mnt/chromeos/MyFiles/Downloads B=`basename $0 .sh`
Adr=("-1900 0" "0 1999")
AdrInc=(100 100)
#Maybe thicker lines every 800 numbers.... but wouldn't make sense on tilted axes anyway...
#Anyway, my goal isn't to draw maps, but specify a projection... that will need an affine transformation to account
Axis=(X Y)
axis=(x y)
for i in 0 1
do
    set -- +proj=omerc \
	+lonc=${O[0]} +lat_0=${O[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
for i in 0 1
do
    perl wkt.pl ${Adr[@]} ${AdrInc[$i]} $i|
	gmt mapproject -I -J"${projs[$i]}"
done |
    #https://www.openstreetmap.org/relation/123423
    # https://gist.github.com/Michael-fore/b40d2c0333c2bb2fe715f6b3adf9721c https://www.youtube.com/watch?v=fRTHshCj-L0
    #[out:json][timeout:25]; relation(123423);(._;>;);out body;
    #[out:json][timeout:25]; relation["wikipedia"="en:Winnetka, Illinois"];(._;>;);out body;
    gmt spatial -T$bdyfile |
    gmt 2kml -K -Fl -Wthick,red
echo "${O[@]} (${offset_numbers[@]})"| gmt 2kml -K -O -Nt
gmt 2kml -O -Fl -W5p,blue $bdyfile #gpxsee shows color, viking doesn't, but can be run via shell.
exit
{
echo $place house addressing grid definition estimation by Dan Jacobson.
echo See also village Code § ...
if test "${projs[0]}" = "${projs[1]}"
then gdalsrsinfo "${projs[0]}"
else
    for i in 0 1
do
    echo === Axis ${axis[$i]} ====
    gdalsrsinfo "${projs[$i]}"
    done
    fi
}|perl -pwle 's/$/\r/'|tee glnAdrGrd.txt
exit
# https://stackoverflow.com/questions/60856714/how-to-use-overpass-api