### Makefile
N=Northbrook IL Address Grid version 2
DESC1=Using PLSS data.
## Author: Dan Jacobson https://www.jidanni.org/
## Copyright: https://www.gnu.org/licenses/gpl.html
## Created: 2024-04-26T01:47:38+0000
## Last-Updated: 2024-05-06T22:42:43+0000
##     Update #: 300
##
nb_c.vik:
U = ../../../../../utilities
include $U/m1.makefile
export PERLLIB = $U
export PATH := $U:$(PATH)
L=~/Downloads
B=$L/PLSSFirstDivision.geojsonl
T=northbrook_il
twsp=IL030420N0120E0#Northfield Township
twsp0=IL030010N0010E0#IL03 origin
D=~/IL_CadNSDI_V2.gdb.zip
## Below we just used section centers. So never mind

## The following procedure, not eventually used in our "final
## product", is for when dealing with just a sheet of one mile dots:
## First make a bounding box that contains only the nice dots for our area:
nn=42.151836 #new north, below a line with bad dots
ne=-87.768054 $(nn) #42.156895
sw=-87.898700 42.070233
corners.vik:
corners.123.vik: #check for hidden data
corners.123.csv:corners.csv
	sort -t , -k 1n,1 -k 2n,2 $<|\
	perl -F, -anwle '$$F[-1]=$$.; print join ",", @F' > $@
corners.csv:$D
	ogrinfo -so $< PLSSPoint | grep Column
	ogr2ogr -sql "SELECT SHAPE,OBJECTID FROM PLSSPoint" \
	-spat $(sw) $(ne) \
	-lco GEOMETRY=AS_XY -f CSV /vsistdout/ $< | sed 1d | sort -t , -k 3,3 -o $@
## Now we by hand looking at viking, identify where (FIDs of the) three dots are in our township
fid_dot.csv:
	for i in 56348,700200 26807,100200 26803,100600; \
	do echo $$i|sed s/,/,$(twsp)_/; done|sort -t , -k 1 -o $@
seed.pid.csv:fid_dot.csv corners.csv
	join -t , -1 1 -2 3 -o 2.1 2.2 1.2 $^ > $@
seed.miles.rough.txt:seed.pid.csv
	perl -MTwsp2mi -we Twsp2mi::pid2gcp_mi $< > $@
corners.miles.rough.csv:corners.csv seed.miles.rough.txt
	set -x; perl -F, -anwle 'print qq(@F[0,1])' $<|\
	gdaltransform -i -output_xy $$(cat seed.miles.rough.txt)|\
	tr ' ' , > $@
corners.miles.smooth.csv:corners.miles.rough.csv
	perl -F, -anwle '$$_=sprintf "%.0f", $$_ for @F; print join ",", @F;' $< > $@
corners.pid.smooth.csv:corners.miles.smooth.csv
	perl -MTwsp2mi -F, -anwle \
	'print Twsp2mi::point_id_plus_miles(q($(twsp0)_100100),@F)' $< > $@
Grid=origin =>{id => "IL030420N0120E0_300520",\
     address => [ -2800, -1400 ]},num_per_mile => [ (800) x 2 ]
corners.adr.smooth.csv:corners.pid.smooth.csv
	perl -MPointId2Address -nlwe \
	'@K=PointId2Address::id2addr_raw$(\
	)({$(Grid),target=>{id=>$$_}}); print join ",", @K' $< > $@
corners.adr.smooth.check:corners.adr.smooth.csv
	sort $< | uniq --repeated #Well I can choose to ignore it...
gcps_2.txt:corners.adr.smooth.csv corners.csv
	paste -d , $^|perl -F, -anwle 'pop @F; print qq(-gcp @F);' > $@
#OK, now we have all the GCPs ready!

#Now let's make some targets. Note: in this city we want 0 to be 0S and 0W, not 0E and 0N!
EW=( abs $$x ) . ( $$x <= 0 ? "W" : "E" )
NS=( abs $$y ) . ( $$y <= 0 ? "S" : "N" )
i=400
nb.adr.targets: #rows of x,y,label...
	perl -wle \
	'for($$x=-4400; $$x<=0; $$x+=$i){print join ",",   $$x, -50, $(EW)}$(\
	)for($$y=-2000; $$y<=0; $$y+=$i){print join ",", -4350, $$y, $(NS)}' > $@
nb_2.vik:
clean:; rm *.csv *.kmz *.txt *.targets

###### Now just try using section centers straight from the BLM
include $U/m2.makefile
$L/$(twsp).json:
	perl -wle 'for(1..24){push @o, sprintf "$(twsp)SN%02d0", $$_;}$(\
	)print join "%7C", @o;'|fetch_blm_sections > $@
centers.csv:$L/$(twsp).json
	perl -MPLSS_corners -wle '%h=PLSS_corners::centers; for(keys %h){my @o;push @o,$(\
	)@{$$h{$$_}}{qw/lon lat/},$$_ ; print join ",", @o;}' $< > $@
ids.csv:centers.csv
	perl -MBou2colrow -F, -anwle \
	'print join ",", @F[0,1], Bou2colrow::trsec_dir2id($$F[-1],"C");' $< > $@
gcps_c.txt:ids.csv
	perl -F, -MPointId2Address -anlwe \
	'@K=PointId2Address::id2addr_raw$(\
	)({$(Grid),target=>{id=>$$F[-1]}}); print qq(-gcp @K @F[0,1])' $< > $@
nb%.lonlat.csv:nb.adr.targets gcps%.txt
	perl -F, -anwle 'print qq(@F[0,1])' $< | gdaltransform -output_xy \
	$$(cat gcps$*.txt) | tr \  , > $@
nb%.csv:nb%.lonlat.csv nb.adr.targets
	paste -d , $^|perl -F, -anwle 'print join ",",@F[0,1,-1]' > $@
nb_c.vik:
diff:nb_c.csv nb_2.csv
	paste -d , $^ | perl -F, -awnle 'print qq(@F[1,0,4,3])' |\
        geod +ellps=WGS84 -I|sort -k 3n,3 #1 to 6 meter differences. OK!