#!/bin/sh
set -e

# Before running this script you need to geocode member locations, for
# example by simply joining a list of member zips with zip
# coordinates:

#$ curl http://www.erikbolstad.no/nedlasting/postnummer.txt > postnummer.txt
#$ cat postnummer.txt | cut -f 1,7,8|tr "\t" "|" > postnummer.list
#$ perl -ne 'print if /^\d\d\d\d\n/' memberzips.txt > memberzips.filtered.list
#$ sqlite3 memberzips.db "create table members (zip text); create table zips (zip text, lat float, lon float);"
#$ sqlite3 memberzips.db ".import memberzips.filtered.list members"
#$ sqlite3 memberzips.db ".import postnummer.list zips"
#$ sqlite3 memberzips.db "select zip,lon,lat from members join zips using (zip)" > memberlocations.txt

# Then import them into an EPSG:4326 grass location called nuugmembers:
#> v.in.ascii input="memberlocations.txt" output="memberlocations" format="point" fs="|" skip=0 columns="zip varchar(4), x double precision, y double precision" x=2 y=3 z=0 cat=0

# Download r.out.gmap:
#$ wget http://svn.osgeo.org/grass/grass-addons/raster/r.out.gmap/r.out.gmap
#$ chmod +x r.out.gmap

# Then run this script in batch mode:
#$ GRASS_BATCH_JOB=/path/to/this/script /path/to/bin/grass -text /path/to/grassdata/norway_ll/nuugmembers

# The files will be put in output/. To rename them to the structure
# used by openlayers use the rename command bundled with perl:
#$ (cd output; rename -f 's!z(\d+)_(\d+)x(\d+).png!$1/$2/$3.png!; mkdir "$1" unless -d "$1"; mkdir "$1/$2" unless -d "$1/$2"' *.png )

### End of documentation

# Create google-style spherical mercator location and a mapset
g.proj -c location=norway_merc proj4="+proj=merc +lat_ts=0.0000000000 +lon_0=0.0000000000 +k_0=1.0000000000 +x_0=0 +y_0=0 +no_defs +a=6378137 +b=6378137 +nadgrids=@null"
g.mapset mapset=PERMANENT location=norway_merc
g.region n=11543712 s=7940625 e=3540851 w=446986 res=256 -s
g.mapset -c mapset=nuugmembers
g.region -d

# Reproject member coordinates
v.proj input="memberlocations" location="norway_ll" mapset="nuugmembers"

# Create temp and output dirs
mkdir -p color alpha output

for z in 5 6 7 8 9; do
    if [ $z = 5 ]; then
	stddeviation=15000;
	mult=100000
    elif [ $z = 6 ]; then
	stddeviation=15000
	mult=75000
    elif [ $z = 7 ]; then
	stddeviation=10000
	mult=50000
    elif [ $z = 8 ]; then
	stddeviation=5000
	mult=25000
    elif [ $z = 9 ]; then
	stddeviation=2500
	mult=12500
    fi

    g.region -d

    # Set optimum resolution for zoom level (code stolen from r.out.gmap)
    PX2M=$(echo "6378137 / 5340353.715440871795" | bc -l) #correction from meters to pixels(zoom 17)
    res=$(echo "(2^(17-$z)) * $PX2M" | bc -l)
    g.region res=$res

    # Generate density map
    v.kernel input=memberlocations output=memberdensity_z$z stddeviation=$stddeviation mult=$mult
    r.mapcalc "calibrated_z$z=int(graph(memberdensity_z$z, 0,0, 0.000001,1, 0.00005,50, 0.00025,100, 0.001,150, 0.01,200, 0.014,255, 255,255))"

    #r.colors calibrated color=rules <<EOF
#0 white
#50 245:245:30
#150 100:230:25
#175 100:230:25
#255 black
#end
#EOF
    #r.colors calibrated color=rules <<EOF
#0 white
#50 green
#255 black
#end
#EOF
    r.colors calibrated_z$z color=rules <<EOF
0 white
50 15:0:235
150 orange
255 black
EOF

    # Generate opacity mask
    g.copy rast=calibrated_z$z,alpha_z$z --o

    r.colors alpha_z$z color=rules <<EOF
0 white
25 50:50:50
255 50:50:50
end
EOF

    # Output tiles
    ./r.out.gmap input=calibrated_z$z zoom=$z interp=cubic outdir=color
    ./r.out.gmap input=alpha_z$z zoom=$z interp=cubic outdir=alpha

    # Add opacity mask
    for file in `ls color/*.png`; do
	tile=`basename $file`;
	convert color/$tile -channel RGB -separate alpha/$tile -channel RGBA -combine output/$tile;
    done
done