DS=-dsco NAME="Youshi Industrial Park grid addressing and street naming exercise" \
-dsco DESCRIPTION="See $(subst $(HOME)/,https://www.,$(PWD))"
## Author: Dan Jacobson https://www.jidanni.org/
## Copyright: https://www.gnu.org/licenses/gpl.html
## Created: 2026-04-03T22:00:51+0000
## Last-Updated: 2026-04-09T11:26:36+0000
##     Update #: 357
d=~/Downloads/d/data
road84=$d/youshi_roads.geojson
addr84=$d/youshi_addrs.geojson
look:
include ../../../../utilities/integer_graticule.makefile
MAPS= graticule_k.gpx houses.gpx #gcps.gpx #vlr.gpx #points.gpx
MAPS= houses.gpx graticule_k.gpx
m0=cross_25_75.gpx house_2437_2437_7459.gpx house_7454_1756_7455.gpx
MAPS=graticule_k.kmz single_axis.gpx
youshi.kmz:graticule_k.kml single_axis.kml $(m0:gpx=kml)
	ogrmerge -f LIBKML -overwrite_ds -o $@ $^ $(DS)
_x = x0=1400 xi=100 x1=3400
_y = y0=6900 yi=100 y1=7700
#Blow up:
#_x = x0=1750 xi=5 x1=1800
#_y = y0=7420 yi=5 y1=7477

wt=WHERE \"addr:housenumber\" LIKE '%號'
O=ogrmerge -overwrite_ds -o $@ $^
addr.gpkg:$(addr84) reverse_gcps.txt
	$R --optfile reverse_gcps.txt -tps $S \
	"SELECT DISTINCT \"addr:street\","\
	"\"addr:housenumber\", geometry AS geom FROM "\
	"$(notdir $N) WHERE \"addr:housenumber\" LIKE '%號'"
original_addresses.gpkg:addr.gpkg
	$R -nlt POINT $S "SELECT DISTINCT printf('%s%c%s',\"addr:street\","\
	"CHAR(10), \"addr:housenumber\") AS name, geom FROM $(notdir $N) $(wt)"
original_addresses84.gpx:$(addr84)
	$R -nlt POINT $S "SELECT DISTINCT printf('%s%c%s',\"addr:street\","\
	"CHAR(10), \"addr:housenumber\") AS name, geometry FROM $(notdir $N) $(wt)"
houses.gpkg:addr.gpkg
	$R $S "SELECT DISTINCT printf('%d-%d%c%d|%d',floor(x(geom)),ceiling(x(geom)),"\
	"CHAR(10),floor(y(geom)),ceiling(y(geom))) AS name, geom FROM $(notdir $N)" -nlt POINT
gcps.gpx:gcps.txt
	$R -if CSV -nlt POINT -oo X_POSSIBLE_NAMES=field_4 -oo Y_POSSIBLE_NAMES=field_5 $S \
	"SELECT geometry, field_2 || '/' || field_3 AS name FROM $N"
gcps.gpkg:gcps.txt
	$R -if CSV -nlt POINT -oo X_POSSIBLE_NAMES=field_4 -oo Y_POSSIBLE_NAMES=field_5 $S \
	"SELECT geometry AS geom, field_2 || '/' || field_3 AS name FROM $N"
gcps_local.gpkg:gcps.txt
	$R -if CSV -nlt POINT -oo X_POSSIBLE_NAMES=field_2 -oo Y_POSSIBLE_NAMES=field_3 $S \
	"SELECT geometry AS geom, FORMAT('%.0f/%.0f', field_2/100, field_3/100) AS name FROM $N"
points.gpkg:$(road84)
	$R -explodecollections -nlt POINTS $S \
	"SELECT ST_DissolvePoints(GEOMETRY) AS geom,name FROM $(notdir $N)"
pv.vrt:points.gpkg road_values.csv; $O
vlr.gpkg:pv.vrt
	$R $S "SELECT road_values.value AS name, geom FROM road_values,points "\
	"WHERE road_values.name = points.name ORDER BY geom"
ints.gpkg:vlr.gpkg
	$R $S "SELECT t1.geom as geom, t1.name AS n1, t2.name AS n2 "\
	"FROM $N AS t1, $N AS t2 WHERE t1.geom=t2.geom AND NOT t1.name = t2.name GROUP BY t1.geom"
reverse_gcps.txt:gcps.txt
	ogr2ogr -if CSV -of CSV /vsistdout/ $< \
	-lco SEPARATOR=SPACE -lco STRING_QUOTING=IF_NEEDED -sql \
	"SELECT field_1, field_4, field_5, field_2, field_3 FROM $N" | grep -- -gcp > $@
gcps.txt:ints.gpkg
	ogr2ogr -of CSV /vsistdout/ $< -lco STRING_QUOTING=IF_NEEDED -lco SEPARATOR=SPACE $S \
	"SELECT '-gcp', min(n1, n2) AS x, max(n1, n2) AS y,"\
	"x(geom) AS lon, y(geom) AS lat FROM $N"|\
	perl -nwle 'print if /\d/' > $@
%.gpx:%.gpkg gcps.txt
	set -eu; set -- $$(ogrinfo -q $<); shift 2; case $$@ in \(Point\)) n=POINT;;\
	*Multi\ Line\ String*) n=MULTILINESTRING;; *Line*) n=LINESTRING;;\
	*Polygon*|*Multi*) echo No $$@ allowed in GPX! 1>&2; exit 67;; \
	*) echo $< --\> $@: Not ready for type \"$$@\" 1>&2; exit 54;; esac;\
	$R -nlt $$n $S 'SELECT name, geom AS geometry FROM "$N"' --optfile gcps.txt -tps
%.kml:%.gpkg gcps.txt
	$R $S 'SELECT name, geom AS geometry FROM "$N"' \
	-of LIBKML --optfile gcps.txt -tps
look:$(MAPS); > /tmp/$(VIEWER)errs.txt $(VIEWER) $^ & sleep 11h
S=-dialect SQLite -sql
R=ogr2ogr $@ $< -nln $(basename $@)
N=$(basename $<)
VIEWER=viking
MAKEFLAGS+=--warn-undefined-variables -r #-B
house_%.gpkg: #Usage: house_Label_X_Y.gpkg
	set -ue -- `echo $*|tr _ ' '`; test $$# -eq 3; \
	ogr2ogr -nlt POINT -nln $(basename $@) $@ :memory: $S \
	"SELECT ST_Point($$2,$$3) AS geom, $$1 AS name"
cross_25_75.gpkg:graticule_k.gpkg #to make the picture of the false origin intersection
	$R $S "SELECT * FROM $N WHERE name % 2500 = 0"
clean:; rm -vf *.csv *.txt *.gpkg *.gpx *.vrt *.kml
names.csv:$(road84) #For our name list below.
	$R $S "SELECT DISTINCT name FROM $(notdir $N) WHERE name like '%路' ORDER BY name"
#We are only doing the main part of the industrial area, not roads called 東
#Overpass-Turbo queries used:
# way["highway"]["name"~"^(幼|工|青)"]({{bbox}});
# node["addr:street"~"^(幼|工|青)"]({{bbox}});
road_values.csv:
	set -ue -- $(road_values); { echo name,value; while test $$# -gt 0; do \
	case $$2 in [0-9]*[0-9]) echo $$1,$$2;; esac; shift 2; done;} > $@
#x: would distort grid:
road_values=\
青一路 1400\
幼一路 1600\
青二路 1600\
青三路 1750\
幼二路 1900\
幼三路 1950\
幼四路 2200\
幼五路 2500\
幼六路 2800\
幼七路 3100\
幼九路 3400\
工十路 6900\
工六路 x7100\
青一路3巷 7100\
青二路3巷 7100\
青三路3巷 7100\
青年路 7200\
工四路 7300\
青一路8巷 7300\
青二路8巷 7300\
青三路8巷 7300\
工二路 7400\
青一路11巷 7400\
青二路11巷 7400\
青三路11巷 7400\
幼獅路 7500\
工一路 x7600\
工九路 7600\
工三路 7700\
工七路 7700
interval=100
bounding_box:gcps.txt #In case I'm curious...
	ogr2ogr -of CSV -if CSV /vsistdout/ $< $S \
	"SELECT FORMAT('_x = x0=%d xi=%d x1=%d%c_y = y0=%d yi=%d y1=%d',"\
	"min(field_2), $(interval), max(field_2), char(10), min(field_3),"\
	"$(interval), max(field_3)) FROM $N"|tr -d \"
F=5000
single_axis.gpkg: #Need two lines so labels fall on each one
	ogr2ogr $@ :memory: -nln $(basename $@) -nlt LINESTRING $S \
	"SELECT MakeLine(MakePoint($F-$F,$F),MakePoint($F-$F,$F+$F)) AS geom,'0000' AS name UNION "\
	"SELECT MakeLine(MakePoint($F-$F,$F),MakePoint($F,$F)) AS geom, $F AS name UNION "\
	"SELECT MakeLine(MakePoint($F,$F),MakePoint($F,$F+$F)) AS geom, $F AS name"