#!/bin/bash
# usage:
#   geotiff2kal_2 map.tif > map.kal
#
# convert GeoTIFF Calibration to Glopus KAL
#
# frank334, 9.2.2008
#
#
#   test mit bernahme der Grad-Koordinaten

if [ $# -ne 1 ]; then
  echo 'Syntax: geotiff2kal_2 map.tif > map.kal'
  exit
fi

image=$1
imageinfo=$1.info
gdalinfo "$image" > "$imageinfo"

base=${image%.*}
tempcoo=$base.coordinates

#utmzone=`cat "$imageinfo" | grep "^PROJCS" | sed "s/^PROJCS.*UTM zone \(.*\)\".*/\1/"`
datum=`cat "$imageinfo" | grep "GEOGCS" | sed "s/GEOGCS\[\"\(.*\)\".*/\1/"`
datum_epsg=`cat "$imageinfo" | grep "^        AUTHORITY.*,$" | sed "s/.*AUTHORITY\[\"EPSG\",\"\(.*\)\".*/\1/"`

sizex=`cat "$imageinfo" | grep "^Size is " | sed "s/^Size is \(.*\), \(.*\)/\1/"`
sizey=`cat "$imageinfo" | grep "^Size is " | sed "s/^Size is \(.*\), \(.*\)/\2/"`

# extract coordinates
ul=`cat "$imageinfo" | grep "^Upper Left" | sed "s/^Upper Left  (.*) ( \(.*\), \(.*\))$/\1 \2/"`
ll=`cat "$imageinfo" | grep "^Lower Left" | sed "s/^Lower Left  (.*) ( \(.*\), \(.*\))$/\1 \2/"`
ur=`cat "$imageinfo" | grep "^Upper Right" | sed "s/^Upper Right (.*) ( \(.*\), \(.*\))$/\1 \2/"`
lr=`cat "$imageinfo" | grep "^Lower Right" | sed "s/^Lower Right (.*) ( \(.*\), \(.*\))$/\1 \2/"`
c=`cat "$imageinfo" | grep "^Center" | sed "s/^Center      (.*) ( \(.*\), \(.*\))$/\1 \2/"`

# make degree format proj-compatible
# 5d6'6.37"E 55d6'48.53"N
ul=`echo $ul | sed "s/\([. ]*\"[EW]\)[ ]*\(.*[NS]\)/\1,\2/" | tr -d " " | tr "," " "`
ll=`echo $ll | sed "s/\([. ]*\"[EW]\)[ ]*\(.*[NS]\)/\1,\2/" | tr -d " " | tr "," " "`
ur=`echo $ur | sed "s/\([. ]*\"[EW]\)[ ]*\(.*[NS]\)/\1,\2/" | tr -d " " | tr "," " "`
lr=`echo $lr | sed "s/\([. ]*\"[EW]\)[ ]*\(.*[NS]\)/\1,\2/" | tr -d " " | tr "," " "`

# DHDN
#<4314> +proj=longlat +ellps=bessel +datum=potsdam +no_defs  <>

# DHDN / Gauss-Kruger zone 2
#<31466> +proj=tmerc +lat_0=0 +lon_0=6 +k=1.000000 +x_0=2500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs  <>
# DHDN / Gauss-Kruger zone 3
#<31467> +proj=tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs  <>
# DHDN / Gauss-Kruger zone 4
#<31468> +proj=tmerc +lat_0=0 +lon_0=12 +k=1.000000 +x_0=4500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs  <>
# DHDN / Gauss-Kruger zone 5
#<31469> +proj=tmerc +lat_0=0 +lon_0=15 +k=1.000000 +x_0=5500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs  <>
# convert MagicMaps (Gauss-Krger DHDN Meridian 9) = Zone 3


# PROJ EPSG Codes aus
# /usr/local/share/proj/epsg

# Datums-Formate im Geogrid
# Generalkarte bersicht: ED50
# Generalkarte: DHDN

#if [ $datum = "DHDN" ]; then
#  echo "DHDN mode (Zone 3, Meridian 9)"
#  epsgcode=31467
#  epsgcode=$datum_epsg
#  sourcedatum="+proj=tmerc +init=epsg:$epsgcode"
#echo $sourcedatum
#elif [ $datum = "ED50" ]; then
#  echo "ED50 mode"
#  echo "EPSG Datum $datum_epsg"
#  sourcedatum="+proj=lcc +init=epsg:$datum_epsg"
#fi

sourcedatum="+proj=longlat"

rm "$imageinfo"

# cs2cs input
# echo "5d6'6.37\"E 55d6'48.53\"N" | cs2cs +proj=longlat +ellps=bessel +datum=potsdam +to +proj=longlat +datum=WGS84 -f %.9f
# 

#echo -e $ul "\n" $ll "\n" $ur "\n" $lr

echo -e $ul "\n" $ll "\n" $ur "\n" $lr | \
cs2cs $sourcedatum +to +proj=longlat +datum=WGS84 -f %.9f | cut -s -f 1 -d " " > $tempcoo
#+datum=$datum

#cat $tempcoo

echo -e 0 0 "\n" 0 $sizey "\n" $sizex 0 "\n" $sizex $sizey | \
paste $tempcoo - | asc2kal "$image" $sizex $sizey

rm $tempcoo
