[PROJ] Help with transformation of GPS coordinates to TK25

[Dominik Vogt]

I've been trying to write a script (Linux command line) that
translates GPS coordinates coming from an Android phone to TK25
sheets in Germany.  As far as I understand:

 * The Android GPS uses epsg:4326 which is also used by Openstreetmap.
   Well, at least putting the coordinates from the device in a marker
   to osm results in the expected position.

   https://www.openstreetmap.org/?mlat=...&mlon=...

 * There is a German page that can do various transformations (for
   botanical mappibg), and the author says that he manually
   compared the borders of the TK quadrants in the online map with
   his printed sheets on paper.

     https://www.orchids.de/haynold/koordinatenermittler2/

   (First select "Karten usw.", then check "Meßtischblatt" and
   "Quadranten", then click "ausblenden".  This show the grid when
   zooming into the map enough.  Second, press the "K" button on the
   left and select "Meßtischblatt-Quadrant".  The pointer should now
   show the TK25 quadrant it's pointing into.

 * The author also says that the TK quadrants are calculated from
   epsg:4314 coordinates.

 * As far as I understand, TK25 divides each degree into ten sheets
   in W-E-direction and six in N-S-direction and is aligned to whole
   degrees.

Are these assumptions correct?

--

So, I use this shell code (zsh) to do the conversion:

-- snip --
bessel_to_tk () {
	local LAT
	local LON
	local TKLAT
	local TKLON
	local TK
	local Q

	while read LAT LON REST; do
		if echo "$LAT$LON" | grep -q "[^0-9.]"; then
			echo "N/A"
		else
			TKLAT="$[560 - 10 * LAT]"
			TKLON="$[6 * LON - 34]"
			TK="$[${TKLAT%.*} * 100 + ${TKLON%.*}]"
			Q=1
			test x"${TKLAT:#*.[5-9]*}" = x && Q="$[Q + 2]"
			test x"${TKLON:#*.[5-9]*}" = x && Q="$[Q + 1]"
			echo "$TK/$Q"
		fi
	done
}

cs2cs +init=epsg:4326 +to +init=epsg:4314 -f "%.12f" |
bessel_to_tk
-- snip --

This fails near borders; for example the north of the Kissinger
Heide near Augsburg is in TK25 7731/2 according to the mentioned
web page by more than 250 m:

  https://www.openstreetmap.org/?mlat=48.2958411471918&mlon=10.9570727683604

But the conversion puts it in 7631/4:

  https://www.openstreetmap.org/?mlat=48.2958411471918&mlon=10.9570727683604

 $ echo 48.2958411471918 10.9570727683604 | ./script
 7631/4

So, _where_ is the offset coming from?  I really don't know how to
fix this?

--

 $ proj
 Rel. 7.2.1, January 1st, 2021

Ciao

Dominik ^_^  ^_^

--

Dominik Vogt