Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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')
share|improve this question
    
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
1  
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
1  
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
show 3 more comments

1 Answer 1

up vote 1 down vote accepted

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.