Add Raster readme.md

requirements

• gdal (and python bindings)
• cartodb up and running
  • it would be nice to use ghost tables branch (CDB-2572) so it can get tables into the dashboard. requires resque to be running also.

quick start

# Download shaded relief from natural earth (http://www.naturalearthdata.com/downloads/10m-raster-data/10m-shaded-relief/)
make prepare
# add the table and overviews
./add_raster.sh example/SR_HR.tif `DBUSER` `DBNAME`

Some notes

Raster related views

${USER} and publicuser requires read permission on raster related views:

GRANT SELECT ON TABLE public.raster_overviews TO publicuser;
GRANT SELECT ON TABLE public.raster_columns TO publicuser;

add_raster.sh

#!/bin/sh

if [[ $# -ne 3 ]]
then
    echo "usage: ${0} <filename.tif> <dbuser> <dbname>"
    exit 1
fi

function scale() {
    if [[ $# -ne 1 ]]
    then
        echo "usage: scale webmercator_filename.tif"
        exit 1
    fi
    local _FILE=$1
    local input_scale=`gdalinfo ${_FILE} | grep '^Pixel Size' | sed 's/.*(\([^,]*\).*/\1/'`
    local z0=156543.03515625
    local factor=`echo "scale=10; $z0/$input_scale" | bc`
    local pow2=`echo "pw=l($factor)/l(2); scale=0; pw/1" | bc -l`
    local out_scale=`echo "scale=10; $z0/(2^$pow2)"| bc`
    local out_scale=`echo "scale=10; $out_scale-($out_scale*0.0001)" | bc`
    echo ${out_scale}
}

ORIG_FILEPATH=$1
DBUSER=$2  # "development_cartodb_user_4293f3c5-b2e3-4d70-9ff3-61ca4ad5747f"
DATABASE=$3  # "cartodb_dev_user_4293f3c5-b2e3-4d70-9ff3-61ca4ad5747f_db"


FILENAME=`basename $1`
_FILENAME="${FILENAME%.*}"
WEBMERCATOR_FILENAME="${_FILENAME}_webmercator.tif"
ALIGNED_WEBMERCATOR_FILENAME="${_FILENAME}_aligned_webmercator.tif"
TABLE_NAME="${_FILENAME}"


gdalwarp -t_srs EPSG:3857 ${ORIG_FILEPATH} tmp/${WEBMERCATOR_FILENAME}

OVERVIEWS=`python raster_overviews.py ${ORIG_FILEPATH}`

OUTPUT_SCALE=`scale tmp/${WEBMERCATOR_FILENAME}`

gdalwarp -tr ${OUTPUT_SCALE} -${OUTPUT_SCALE} tmp/${WEBMERCATOR_FILENAME} tmp/${ALIGNED_WEBMERCATOR_FILENAME}

raster2pgsql -s 3857 -t 128x128 -C -Y -I -f the_raster_webmercator -l ${OVERVIEWS} tmp/${ALIGNED_WEBMERCATOR_FILENAME} ${TABLE_NAME} > out/${TABLE_NAME}.sql


psql -U${DBUSER} -d ${DATABASE} -f out/${TABLE_NAME}.sql
psql -U${DBUSER} -d ${DATABASE} -c "select cdb_cartodbfytable('${TABLE_NAME}');"

Makefile

prepare:
	chmod u+x add_raster.sh
	mkdir -p example
	mkdir -p out
	mkdir -p tmp
	# curl -L http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/raster/SR_HR.zip -o example/SR_HR.zip
	cd example; unzip SR_HR.zip; rm *.html *.txt *.tfw *.zip *.prj

raster_overviews.py

import gdal
import math
import sys

if len(sys.argv) < 2:
	print "usage"
	exit(1)

filename = sys.argv[1]
src = gdal.Open(filename)
bigger_size = max(src.RasterXSize, src.RasterYSize)
max_power = int(math.ceil(math.log(bigger_size/256,2)))
overviews = map(lambda x: int(math.pow(2, x)), range(1, max_power + 1))
overviews = filter(lambda x: x <= 1000, overviews)

print ",".join(map(str, overviews))

overviews = filter(lambda x: x <= 1000, overviews)
that 1000 is because the tool doesn't allows values bigger than 1k