#!/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/" | tail -1`

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

# TOP25:
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/"`


# 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-Krüger 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"
elif [ $datum = "ED50" ]; then
  #echo "ED50 mode"
  #echo "EPSG Datum $datum_epsg"
  sourcedatum="+proj=lcc +init=epsg:$datum_epsg"
fi

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