#!/bin/bash
# Glencoe Illinois USA address grids --- second attempt,
# using GMT's built in grids...
# Making something like
# https://mapwarper.net/maps/75719#Preview_Rectified_Map_tab
# "Let's quit while we're ahead."
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.html
# Created: 2023-07-08T20:45:39+0000
# Last-Updated: 2024-02-13T13:09:35+0000
#     Update #: 334
set -ue
place=Glencoe\ IL\ USA
bdyfile=glencoeBdy.gmt
if ! test -s $bdyfile
then
    ogr2ogr -f GMT -where "OGR_GEOMETRY='POLYGON'" $bdyfile glencoeBdy.kml
fi
#We do X and Y independently.
#They might need different +k_0, and they might not be at 90 angles in some cities...
set -eu
#set ESRI:102445 # EPSG:3078 #gdalsrsinfo "$*"
#gdalsrsinfo "$*"
# https://codelibrary.amlegal.com/codes/glencoe/latest/glencoe_il/0-0-0-19602
# § 30-65 BASE LINES DESIGNATED; STARTING NUMBERS; DIRECTIONS OF NUMERICAL INCREASE.
# We note that they are talking about what is now OLD Green Bay Rd.!

L=(-87.7415807 42.1192853) # https://communitymapviewer.gisconsortium.org/GlencoeIL
# and at right angles to the right-of-way of the Chicago and
# Northwestern Railway Company
#: Here we get two points on the west side of that ROW:
RRs=(-87.7447873 42.1192976)
#: and we can use where it hits Lake Cook Rd... for the north point.
RRn=(-87.7725841 42.1523955)
# : Let's find out where that perpendicular line they mentioned hits the ROW west side
O=($(echo ${L[@]} | gmt project -Frs -C${RRs[0]}/${RRs[1]} -E${RRn[0]}/${RRn[1]})) #our "origin"...
# https://forum.generic-mapping-tools.org/t/how-to-make-mapproject-give-same-azimuth-as-proj-orgs-geod/4125/2
# Get the azimuth of the luckily rather straight line RR ROW
set -- $(echo ${RRs[@]} ${RRn[@]} | gmt mapproject -AF+v -fg -o4)
azimuth=($@ $@) #both the same for this town
num_per_mile=$(gmt math -Q 300 3000 5280 DIV DIV =) #num / ft /  (ft/mi) = num / mi
#If I want +k_0 applied to +x_0, +y_0, then I must apply it myself beforehand.
#
# 300 nums...
#x_0=$(gmt math -Q -300 10 MUL 5280 DIV $num_per_mile MUL 1609.344 MUL =)
#offset=($(gmt math -Q -3000 5280 DIV $num_per_mile MUL 1609.344 MUL =) 0)
offset=(+x_0=$(gmt math -Q -3000 5280 DIV $num_per_mile MUL 1609.344 MUL =) +y_0=0)
#Let's break that down:
#We need meters of offset.
#But "+units=us-mi" will act upon this, so we need miles of offset.

# -3000 5280 DIV: backwards 3000 ft/5280 (ft/mi) unit: miles
# acted upon it, which multiplies it by 
#  Remember that "+units=us-mi" will affect this, just like if we did
#  1609.344 DIV so we do
#  1609.344 MUL: now.

#then we will have miles
#  $num_per_mile MUL: times house numbers per mile, so now we have house numbers

# spell check all files...
#use our magic formula to adjust Y some meters...

#I bet the 1799 etc. Green Bay Road are from the Winnetka system.

#gmt 2kml man page: make -O docstring match -K, also for usage
D=/mnt/chromeos/MyFiles/Downloads B=`basename $0 .sh`
Adr=("-1500 0" "0 1799")
AdrInc=(100 100)
scale=($num_per_mile $(gmt math -Q $num_per_mile .97 MUL =)) #tweak for Y, which makes it match the ground better for some reason.
#set -x
#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
#for 97 x 100 stuff...
axis=(X Y)
for i in 0 1
do
    set -- +proj=omerc +lonc=${O[0]} +lat_0=${O[1]} +alpha=${azimuth[$i]} +gamma=0 \
	${offset[@]} +units=us-mi +k_0=${scale[$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;
    gmt spatial -T$bdyfile |
    gmt 2kml -K -Fl -Wthick,red >  $D/$B.kml
gmt 2kml -O -Fl -W5p,DarkOrange $bdyfile >> $D/$B.kml
#for i in L RRs RRn O; do eval echo \${$i[@]} $i; done|
#    gmt 2kml -O -Nt >> $D/$B.kml
{
echo $place house addressing grid definition estimation by Dan Jacobson.
echo See also village Code § 30-65
for i in 0 1
do
    echo === Axis ${axis[$i]} ====
    gdalsrsinfo "${projs[$i]}"
done
}|perl -pwle 's/$/\r/'|tee glnAdrGrd.txt