Get highest peaks in each section of Taipower Coordinate System (8000m x 5000m) #mapstew #taipower

highest_peaks_in_taipower_sections.sh

#! /bin/bash

#===============
# Constants
#===============

TARGET_OSM_FILE=$1
OUTPUT_FILE=$2

SPACE_X=8000
SPACE_Y=5000

declare -A sections
sections[A]=170000,2750000
sections[B]=250000,2750000
sections[C]=330000,2750000
sections[D]=170000,2700000
sections[E]=250000,2700000
sections[F]=330000,2700000
sections[G]=170000,2650000
sections[H]=250000,2650000
sections[J]=90000,2600000
sections[K]=170000,2600000
sections[L]=250000,2600000
sections[M]=90000,2550000
sections[N]=170000,2550000
sections[O]=250000,2550000
sections[P]=90000,2500000
sections[Q]=170000,2500000
sections[R]=250000,2500000
sections[T]=170000,2450000
sections[U]=250000,2450000
sections[V]=170000,2400000
sections[W]=250000,2400000

#===============
# Functions
#===============

twd67_to_wgs84() {
  cs2cs -f %.7f \
    +proj=tmerc \
    +lat_0=0 +lon_0=121 \
    +k=0.9999 \
    +x_0=250000 +y_0=0 \
    +ellps=aust_SA \
    +towgs84=-750.739,-359.515,-180.510,0.00003863,0.00001721,0.00000197,0.99998180 \
    +units=m \
    +to \
    +init=EPSG:4326 | \
  tr ' ' '\t'
}

get_highest_peak() {
  echo Filter peaks in boundary: $1 $2 $3 $4 $5
  osmconvert $TARGET_OSM_FILE -b=$2,$3,$4,$5 --csv="name ele" | \
  sort -gk2 -r | head -1 | cut -f1 | \
  sed "s/^/$code,/" >>$OUTPUT_FILE
}

#=============================
# For loop for each boundaries
#=============================

echo code,name >$OUTPUT_FILE

for section in ${!sections[@]}; do
  for xi in {0..9}; do
    for yi in {0..9}; do
      IFS=',' read x y <<<"${sections[$section]}"
      let left=$x+$SPACE_X*$xi
      let bottom=$y+$SPACE_Y*$yi
      let right=$left+$SPACE_X
      let top=$bottom+$SPACE_Y
      code=${section}${xi}X${yi}X

      echo $left $bottom 0 $code
      echo $right $top 0
    done
  done
done | \
twd67_to_wgs84 | \
while { read west south ele code; read east north ele; }; do
  get_highest_peak $code $west $south $east $north
done | \
nl

Makefile

all: peaks.db

taiwan-latest.osm.pbf:
	curl -O http://download.geofabrik.de/asia/$@

TEMP_FILE := temp.o5m
peaks.o5m: taiwan-latest.osm.pbf
	osmconvert $< --drop-ways --drop-relations -o=$(TEMP_FILE)
	osmfilter $(TEMP_FILE) --keep=natural=peak -o=$@
	rm $(TEMP_FILE)

peaks.csv: peaks.o5m
	./highest_peaks_in_taipower_sections.sh $< $@

peaks.db: peaks.csv
	sqlite3 $@ -cmd '.mode csv' ".import $< peaks"