# 曲折單軸門牌號碼規劃 Twisted single axis address planning Makefile
# Project URL     : http://jidanni.org/geo/house_numbering/twist.html
# Author          : 積丹尼 Dan Jacobson http://jidanni.org/
# Copyright       : http://www.fsf.org/copyleft/gpl.html
# Created On      : April 2005
# Last Modified On: Thu Jul 14 04:50:48 2005
# Update Count    : 888

# Main software used:
# http://packages.debian.org/make
# http://packages.debian.org/grass
base_region:#TWD97=WGS84 datum, from orthophotomap.
	g.region n=2680500 s=2678300 w=233000 e=234300 nsres=1 ewres=1
# 台中縣東勢鎮上城社區
# Shangcheng Community, Dongshi Township, Taizhong County, Taiwan

#- 將背景圖掃描後，裝入及矯正
#- Get background orthophotomap into the system and rectified

# After scanning the paper orthophotomap, we do
image_in:
	r.in.gdal -o input=orthophoto.pgm output=rawortho
	i.group group=s2 input=rawortho
	i.target -c group=s2
image_rectify:
	i.points #by hand
	i.rectify -c input=rawortho group=s2 order=1 extension=1
	sleep 5 #as it is mailing in the background
	g.rename rast=rawortho1,orthophoto
	g.remove rast=rawortho

#- 建立等門牌線勢輔助曲線。 等門牌線將與之平行。
#- Make contour trend guides: curves that address contours will run
#- parallel with.

# Create a map with several equal address value contours.
# Just curves, no special categories or values for now. Can add more
# later as needed...
trend:base_region
	v.digit -n map=$@ bgcmd="d.rast map=orthophoto"

#- 建立主軸曲線，與勢線垂直。 將沿此線依距離決定社區整體門牌。
#- Perpendicular to the trend curves, establish the axis curve. The
#- entire community's addresses will be determined according to
#- distance along this axis.

# Now make a map with one curve, e.g., along a road, that will be our
# axis. We start at where house numbers will be counted from 0, and
# number according to distance along this axis. Make sure the axis
# crosses all trend curves.
axis_parts:base_region
	v.digit -n map=$@ bgcmd="d.rast map=orthophoto"

#- 輸入彎端群。 其門牌將由接點獨立發展。
#- Enter bent ends. Their numbers will increase or decrease from where
#- they attach to our address contour surface.

# Now make a map with curves along the bent segments that we will
# number independently continuing from where one of their ends is
# attached to our address "surface".
bends:base_region
	v.digit -n map=$@ bgcmd="d.rast map=orthophoto"
#Start and end points are significant! Start: the place where we begin
#to leave the overall system.

#END OF PRECIOUS HAND ENTERED DATA

#Various steps, done one by one:

M=9999#a large number
ribcage:base_region
	v.category input=trend output=trend1 option=del
#(As I entered the road in two pieces, so need to:)
	v.build.polylines input=axis_parts output=axis
# Make a point for the origin of address numbering
	v.to.points -t input=axis output=n3a dmax=$M
	v.extract input=n3a output=origin layer=2 list=1 #Only need origin
	g.remove vect=n3a
# Combine axis and trend curves
	v.patch input=axis,trend1,origin output=b5
	g.remove vect=trend1,origin,axis
# give curves cat numbers
	v.category type=line input=b5 output=b7
	g.remove vect=b5
# break all into separate curves
	v.clean input=b7 output=b8 tool=break
	g.remove vect=b7
	v.category input=b8 output=n9 option=chlayer layer=1,3
	g.remove vect=b8
	-g.remove vect=$@
	v.category input=n9 output=$@ option=add type=line
	g.remove vect=n9
see_ribcage:base_region #inspect at this stage if desired
	d.erase
	d.vect map=ribcage display=shape,cat llayer=3 #Axis=1 vs. ribs
ribcage.ps:base_region #output to file
	echo 'vlines ribcage; color black; end; vpoints ribcage;\
	color black; fcolor black; end;'|tr \; \\n|\
	ps.map scale=1:40000 output=$@.tmp #but ugly box:
	perl -pwe 's/^\d.* B D$$//#Until I can remove it in grass!' $@.tmp>$@
	rm $@.tmp
%.png:%.ps
	convert $? $@
	convert -trim $@ $@
	convert -monochrome $@ $@
# 得排骨狀。 勢線排骨接主軸脊椎。 全獨骨。
# We now have a rib cage of trend curves attached to a backbone of the axis.
# All are separate bones.

iso:base_region
	-echo DROP TABLE roads|db.execute
	echo CREATE TABLE roads \(cat int,oline int,cost double\)|db.execute
	-v.db.connect map=ribcage table=roads
	-v.db.connect map=ribcage layer=3 table=roads
	v.to.db map=ribcage option=cat
	v.to.db map=ribcage option=length type=line col=cost
	v.to.db map=ribcage qlayer=3 column=oline option=query qcolumn=cat
	echo 'UPDATE roads SET cost=0 where oline<>1'|db.execute #all ribs=0
#I.e., each rib is at the same "height" as where it meets the backbone.
#	echo SELECT \* FROM roads|db.select #to inspect
	v.net.iso input=ribcage output=iso1 ccats=1 \
	costs=$$(seq -s, 10 10 2000) afcolumn=cost #band width 10 meters
	v.category input=iso1 output=i1 option=chlayer layer=1,2
	g.remove vect=iso1
	v.category input=i1 output=$@ option=add
	g.remove vect=i1
# We now have put the rib cage into 10 meter bands.
	-echo DROP TABLE bands|db.execute
	echo CREATE TABLE bands \(cat int,contour int\)|db.execute
	-v.db.connect -o map=iso table=bands
	v.to.db map=iso option=cat
	v.category input=iso option=print layer=1,2|perl -nwe \
	'm/(\d+)\|(\d+)/ && \
	printf "UPDATE bands SET contour=%d WHERE cat=%d;\n",\
	$$2*10,$$1'|db.execute #elevations in meters
#	echo SELECT \* FROM bands|db.select #to inspect
# We now have put the bands into the database.

surface:base_region #not clip_region!
	-g.remove rast=$@
#- 製等門牌面。
#- Create a raster surface of address "elevations". Takes minutes.
	v.surf.rst input=iso elev=$@ dmax=$M zcolumn=contour

clip_region: #Just the meat.
	g.region n=2680136 s=2678983 e=233955 w=233227
F=400/1000#addresses per meter. Note: hardwired elsewhere in this
#program though.
contours:clip_region
	-g.remove vect=$@
#- 製等門牌線。
#- Make contours.
# We make one contour line per 10 house numbers, i.e., 25 meters:
	r.contour input=surface output=$@ step=25
	echo ALTER TABLE $@ ADD COLUMN address int|db.execute
	echo SELECT level FROM $@|db.select -c|perl -nwe \
	'printf "UPDATE $@ SET address=%d where level=%d;\n",\
	$F*$$_,$$_ unless int $$_%125'|db.execute
# We are now also ready to make labels every five curves, i.e.,
# 5*25=125 meters.

#- 接每彎端之一腳至等門牌面。
#- Attach bent ends of roads that are to be numbered with only one end
#- rooted in our surface:
bends2:base_region
	-g.remove ve=$@
	v.to.points input=bends output=$@ dmax=$M
	echo ALTER TABLE $@_2 ADD COLUMN contour double|db.execute
	v.what.rast vector=$@ rast=surface col=contour layer=2
	echo ALTER TABLE $@_2 ADD COLUMN address int|db.execute
	echo SELECT \* FROM $@_2|db.select -c|perl -F\\\| -anlwe \
	'if($$F[2]==0){$$elev=$$F[3];if($$F[1]=~/7|2|5/){$$f=-1}else{$$f=1}};\
	printf "UPDATE $@_2 SET address=%d where cat=%d;\n",\
	($$elev + $$F[2]*$$f)*$F, $$F[0]'|db.execute
# Road segments 7,2,5 are decreasing. $$F[2]==0 at attachment points...
#	echo SELECT \* FROM $@_2|db.select #to inspect

labels:
# Labels for main contours:
	v.label map=contours column=address labels=contourslabels size=33 \
	border=black background=white
# Labels for bent ends:
	v.label map=bends2 layer=2 column=address labels=bends2labels size=22 \
	background=white

.INTERMEDIATE: twist.ps
twist.ps: clip_region
	echo 'greyrast orthophoto;labels contourslabels; end;\
	labels bends2labels; end; vlines contours; end'|tr \; \\n|\
	ps.map scale=1:5000 output=$@
qv:twist.ps #quick view
	gv $?
#- 製圖供承辦員。
#- Produce map for staff.
%.jpg:%.ps
	convert -trim -quality 13 $? $@
	xli $@
view:twist.jpg#view final product

# Extra notes:
#
# 其實整個規劃未牽涉到真正的緯度。我們選擇了一條（略彎的）路作縱軸，
# 等門牌線亦彎，與縱軸維持垂直９０度。
#
# Indeed when we created our contours, we chose a (slightly curving)
# north-south road as axis. At no time were actual northings used. At
# all times distances were measured along curves, and then numbers
# were assigned according to these distances.
#
# 在本社區， 在多數的路，我們犧牲了「每公里四百個號碼」之精準，為了達
# 到「平行的路，號碼同步」之更高目標。
#
# For this community, for most roads, we sacrifice exact 400 numbers
# per kilometer accuracy, for the higher goal of parallel roads having
# parallel house numbers.
#
# 圖上顯得略怪， 但對路中人才有意思。 在北美， 過了第十九橫街，所有的
# 路的門牌就達１９００號了。本案中之東西向路不落於整百位置，但我們仍要
# 求「過某路，門牌同步」。
#
# We always must put ourselves with the man in the street, coming up
# with contours that might look funny on paper, but hopefully will
# make sense at the scene. We intend the same east-west cross road to
# signify reaching the same address value, for addresses on all
# north-south roads crossing it. One would want to further give that
# east-west cross road a name relating to its value, e.g., 19th Street
# for 1900, but for our case we are not dealing with such nice
# hundred-numbers, however the subordinate lane number will naturally
# still contain such information.

# Author's website maintenance, etc.:
web:
	cp -avu twist.jpg ribcage.png ~/jidanni.org/geo/house_numbering/images
	cp -avu Makefile \
	~/jidanni.org/geo/house_numbering/programs/twist/makefile.txt
P=/mnt/usb/pqi
backup:
	-mount $P
	mkdir -p $P/shangcheng
	cp -vua Makefile ~/jidanni.org/geo/house_numbering/twist.html \
	$P/shangcheng
	cd /home/jidanni/tmp/grassdata/shangcheng/addressing && \
	cp -avu group vector/axis_parts vector/bends vector/trend $P/shangcheng