# domain for transverse mercator projection, using proj.4

I have to apologize for the naive question -- it's my first encounter with projections and this software package, and I intend to use it for a non-geographical application (hence I work on a unit sphere) -

when using the transverse mercator transformation, I find that it is not defined for latitude values between pi/2 and 3*pi/2, as well as -pi/2 and -3*pi/2 and thus discontinuous - why should that be, and how do I make this compatible with "my definition" of spherical coordinates, as (r,theta,phi) with -pi<=theta<=pi and -pi/2<=phi<=pi/2? It's easy to see that for phi I can just substract pi/2 (thus set the equator to 0), but what about theta?

Thanks a lot, Ben

``````import matplotlib.pyplot as plt
import numpy             as np

from pyproj import Proj

pi=np.pi

tmerc = Proj("+proj=tmerc +a=1")

def f(a,b): return tmerc(a,b,radians=True,errcheck=True)

nice   =[]
naughty=[]

for i in     np.arange(-2*pi,2*pi,.2):
for j in np.arange(-2*pi,2*pi,.2):
try:    a,b=f(i,j) ; nice.append((i,j))
except: naughty.append((i,j))

plt.plot(np.array(nice   )[:,0],np.array(nice   )[:,1],color='green',marker='o',markeredgewidth=0.,linestyle='None')
plt.plot(np.array(naughty)[:,0],np.array(naughty)[:,1],color='red'  ,marker='o',markeredgewidth=0.,linestyle='None')
``````
-
Did you really mean to write "latitude" values between pi/2 and 3pi/2, and not longitude? If longitude was indeed meant, consider that the usual Mercator projection has poles at latitudes +-90 degrees from the Equator, then turn it on its side to get a TM centered on the Prime Meridian: where are the poles now? :-) – whuber Aug 23 '12 at 20:00
thanks, absolutely right -- longitude was meant, and it sort of makes more sense now. but still, how do I project the (continuous) interval [-pi;+pi]x[-pi/2;+pi/2] ?? – bs448c Aug 23 '12 at 20:48
Have you thought of making two maps, one centered at 0 and the other at 180 degrees, and showing them side-by-side? – whuber Aug 23 '12 at 21:21
cheers, so what seems to work is the following: in my function `f(a,b)` I translate all values a<-pi/2 and a>pi/2 to the correct quadrant by substracting or adding pi, respectively, and inverting the sign of the resulting projection coordinates x,y. in this way I get a mapping from [-pi;pi]x[-pi/2;pi/2] to [-pi;pi]x[-pi/2;pi/2] which is continuous and neighbourhood-preserving as long as I take periodicity into account (y_min = y_max). I tested this with calculating the correlation between great circle distance in "sphere space" and euclidean distance in "projection space". – bs448c Aug 24 '12 at 13:43
Yes, that's exactly what I was suggesting. Most GISes just can't handle limiting values of projections properly, so you're lucky such a simple workaround is available. They can't even handle the periodicity at +-pi properly, for heaven's sake :-) -- just check out all the questions we've had here about problems with crossing that meridian. – whuber Aug 24 '12 at 13:47

## 1 Answer

right, so in the end was seems to be the case is the following:

tmerc projects `[-pi;+pi]x[-pi/2;+pi/2]` to a weird non-continuous projection space `ax[-pi/2;+pi/2]` with `a` between `-2pi` and `-3/2*pi` or between `-pi/2` and `+pi/2` or between `+3/2*pi` and `+2pi`. This was verified by the code above.

My goal was to project my continuous range of spherical coordinates to a continuous range of cartesian coordinates which I managed to achieve by translating longitudes (hope I got it right this time) from the "forbidden area" `pi/2 < |theta| < 3/2*pi` to the proper quadrant, by adding or substracting `pi`, respectively; in this case, also the sign of the resulting, projected coordinates had to be inverted. Also, the two maps "have to be displayed next to each other", by adding `pi` to the projected `y`-coordinate in case of the translation.

in python: instead of the function `f(a,b)` above I used something like

``````def f(a,b):

flag=1

if a<-pi/2: a-=pi ; flag=-1
if a> pi/2: a+=pi ; flag=-1

x,y=tmerc(a,b,radians=True,errcheck=True)

x*=flag ; y*=flag
if flag==-1: y+=pi

return x,y-pi/2
``````

If I now compare the pairwise Great Circle Distance of points equally distributed on the surface of a sphere with their (Euclidean) distance in projection space, I could verify that neighbourhood is preserved by the projection (of course there's some distortion). Plus, it's essential that in my case the `y` values are "wrapped around" since there's periodicity at `+-pi`.

Hope that can be of use for someone.

-