#!/bin/sh
# lake_co ---
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.html
# Created: 2023-06-30T19:44:52+0000
# Last-Updated: 2023-07-27T07:27:18+0000
#     Update #: 627
#
#gmt clear all
set -eu
title="Lake Co. IL  house addressing grid"

# https://codelibrary.amlegal.com/codes/lakecountyil/latest/lakecounty_il/0-0-0-8198
# (a)   Standard numbering system. The official address shall follow the
# Standard Street Address System unless the conditions of division
# (A)(1)(b), below, apply. The addressing system is based on a uniform
# grid pattern with its origin beginning at State Street and Madison
# Street in Chicago. The addressing numbers increase as they go west and
# north from that point. The grid number represents the number of miles
# west or north from the base street plus an address number. Lake
# County's beginning reference point is at the southwest corner of the
# county, being 20 miles north of Madison Street and 29 miles west of

# Well more like 19, but...

# State Street in Chicago. Address numbers with a north directional are
# assigned even on the west side and odd on the east side of a street.
# Address numbers with a west directional are assigned even on the north
# side and odd on the south side of a street. Lake County does not
# contain any south or east addresses under this system.

# https://19thcircuitcourt.state.il.us/DocumentCenter/View/433/Map-of-Lake-County-Townships-PDF

# Get control points, at prominent section, township, range corners:
# Made from consulting e.g., openstreetmap, etc. and copying the location of nodes, etc.
# or even coordinates of middles of intersections from Google, etc. satellite Maps.

# Some can even open iD's Measurement Panel, for a temporary node one created.

MMy=42.3277798 MMx=-88.0047754 # M: middle.
NWy=42.4880173 NWx=-88.1993387 # https://www.openstreetmap.org/node/7539142783
NMy=42.4878332 NMx=-88.0031162
NEy=42.4860190 NEx=-87.8252390 # https://www.openstreetmap.org/node/235405212
SWy=42.1542630 SWx=-88.1996240 # https://www.openstreetmap.org/node/261327615
SMy=42.1536527 SMx=-88.0044681 # https://www.openstreetmap.org/node/7879006451
SEy=42.1523558 SEx=-87.7700587 # https://www.openstreetmap.org/node/235545395
MMax=-19000 MMay=32000 #addresses at the center point

# Get the azimuth of the "on the ground grid".
G="geod -f %.7f -p +ellps=GRS80"
set -- $(echo $SMy $SMx $NMy $NMx | $G -I) #no 90 degree shift if I use WM and EM !!
azPole=$1
azEquator=$(gmt math -Q $azPole 90 ADD =) #"oblique equator"
set -- $(seq 0 3|gmt math STDIN 90 MUL $azPole ADD =)
# Our map should extend a little beyond the control points.
# Here we find out the coordinates of how far for each "cardinal direction,"
# E.g., N 13mi from center, E 20mi from center...
set -- $(for AzDist in "$1 13" "$2 15" "$3 13" "$4 13"
	 do
	     $G +units=us-mi <<EOF
    $MMy $MMx $AzDist
EOF
	 done)
n=$1; shift 3; e=$2; shift 3; s=$1; shift 3; w=$2
# A general address calculator, in case I ever need it:
set proj +proj=pipeline \
    +step +proj=omerc +lonc=$MMx +lat_0=$MMy +alpha=$azPole +gamma=0 \
    +step +proj=unitconvert +xy_out=us-mi \
    +step +proj=affine +xoff=$MMax +yoff=$MMay +s11=1000 +s22=1000
P=$@
ff=`basename $0`
#file:///home/jidanni/Downloads/gmt-6.4.0/doc_release/html/pstext.html?highlight=font%20size#description

#perl -nwle 'if(/^..y=(\S+) ...=(\S+)/){print "$2 $1";}' $0 |
#    gmt info -I.1)
#set -- $(perl -nwle 'if(/^..y=(\S+) ...=(\S+)/){print "$2 $1";}' $0 |
 #   gmt info -I.05/.2/.05/.05+R -C)
#rrr=-R$1/$3/$2/$4+r
apm=$(gmt math -Q .0254 12 5280 MUL MUL 1000 DIV =)
PL='printf "%s ag %d\@+000\@+\n", $F[0], abs $F[1]/1000;'
gmt math -T-12000/12500/6000 T = | gmt convert -i+s$apm,+o$MMax |
    perl -nwale "$PL" > /tmp/x.txt
gmt math -T-12000/11700/6000 T = | gmt convert -i+s$apm,+o$MMay |
    perl -nwale "$PL" > /tmp/y.txt
if false
then
for i in x y
do
    perl -e '
    use strict;
    use warnings FATAL => q[all];
    my %info = (
        center => { x => '$MMax', y => '$MMay' },
        x      => [ -12000, 12500, 6000 ],
        y      => [ -12000, 11700, 6000 ],
        NoDir  => 1
    );
    my %dir              = ( x => [qw/E W/], y => [qw/N S/] );
    my $AddessesPerMeter = 1000 / ( .0254 * 12 * 5280 );
    for (
        my $f = @{ $info{'$i'} }[0] ;
        $f <= @{ $info{'$i'} }[1] ;
        $f += @{ $info{'$i'} }[2]
    ) {
        my $g         = $f + $info{center}{'$i'};
        my $direction = $info{NoDir} ? "" : @{ $dir{'$i'} }[ $g < 0 ];
        my $t         = abs $g / 1000;
        my $r         = abs $g % 1000;
        printf "%d ag %d\@+%03d\@+%s\n", $f / $AddessesPerMeter, $t, $r,
          $direction;
    }
				 ' > /tmp/$i.txt
done
fi
if gmt begin $ff png
then
    set -x
    scale=1:800000
    # Soon can use -WM
    rc=-R$w/$s/$e/$n+r
#    rv=$rrr
    rv=$rc
    vv=$(gmt mapproject $rv -Joa$MMx/$MMy/$azEquator/$scale  /dev/null -V -F -C 2>&1|
		 perl -nwle 'if(/Transform/){@F=split; print $F[5]}')
    gmt basemap $rv  -Joa$MMx/$MMy/$azEquator/$scale -Baf -BNE
    gmt coast -W -Na -Glightbrown -Slightblue
    perl -nwle 'if(/^..y=(\S+) ...=(\S+)/){print "$2 $1";}' $0 | gmt plot -N -Si0.2i
    gmt basemap -TdjBR+w0.5i+l
    set -- $(gmt gmtdefaults|perl -anwle '$#F==2 && $F[0]=~/PRIMARY|SECONDARY/ && $F[2]=~s/black/red/ && print "--$F[0]=$F[2]";')
    gmt basemap -R$vv -Jx$scale $@ -Bxc/tmp/x.txt -Byc/tmp/y.txt
    # file:///home/jidanni/Downloads/gmt-6.4.0/doc_release/html/cookbook/options.html#custom-annotations
    echo  0 0 $title | gmt text
    gmt end
    xli $ff.png
fi
# https://github.com/GenericMappingTools/gmt/issues/7541
# file:///home/jidanni/Downloads/gmt-6.4.0/doc_release/html/cookbook/map-projections.html#oblique-mercator-projection-jo-jo