coordinate transformation: Basemap vs. Proj

_Timothy_W_Hilton · May 4, 2010, 8:03pm

Hello matplotlib users,

I am having trouble understanding the coordinate transformations in
Basemap and pyproj. I have gridded MODIS vegetation data, with upper
left corner and lower right corner given in projection coordinates
(meters). I want to contour the data with Basemap. The data are in a
sinusoidal projection, but the coordinates do not correspond to what
Basemap seems to expect.

The code below illustrates the problem. Proj translates the upper
left to lat/lon correctly (-92.327237416031437, 30.141972433747089),
while Basemap does not.

#-------- code --------
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import pyproj

ulm = [-8895604.1573329996, 3335851.5589999999] #upper left, meters
lrm = [-7783653.6376670003, 2223901.0393329998] #lower right, meters

sinu = pyproj.Proj(proj='sinu', lon_0=0.0, x_0=0.0, y_0=0.0)
m = Basemap(projection='sinu', resolution=None, lon_0=0.0)

print "ULM: " + str(ulm)
print "Proj: " + str(sinu(ulm[0], ulm[1], inverse=True))
print "Basemap: " + str(m(ulm[0], ulm[1], inverse=True))
#----- end code --------

This gives:
ULM: [-8895604.1573329996, 3335851.5589999999]
Proj: (-92.327237416031437, 30.141972433747089)
Basemap: (-159.99950210056144, -59.99995206181125)

I'm sure I'm missing something really simple, but I've read a lot of
documentation and I'm not sure what.

Many thanks for any help.

Best,
Tim

···

--

Timothy W. Hilton
PhD Candidate, Department of Meteorology
The Pennsylvania State University
503 Walker Building, University Park, PA 16802
hilton@...3085...

_Jeff_Whitaker1 · May 4, 2010, 10:01pm

Tim: Basemap is using pyproj under the hood, but only supports a subset of possible proj4 projections. The basemap sinusoidal projection is global - you can't specify a subregion of the globe. I think that's where the discrepancy is coming from. I'm sure there's a way to plot your MODIS data on a global sinusoidal projection - but it will involve transforming the coordinates to the Basemap global sinuosidal coordinate system.

-Jeff

···

On 5/4/10 2:03 PM, Timothy W. Hilton wrote:

Hello matplotlib users,

I am having trouble understanding the coordinate transformations in
Basemap and pyproj. I have gridded MODIS vegetation data, with upper
left corner and lower right corner given in projection coordinates
(meters). I want to contour the data with Basemap. The data are in a
sinusoidal projection, but the coordinates do not correspond to what
Basemap seems to expect.

The code below illustrates the problem. Proj translates the upper
left to lat/lon correctly (-92.327237416031437, 30.141972433747089),
while Basemap does not.

#-------- code --------
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import pyproj

ulm = [-8895604.1573329996, 3335851.5589999999] #upper left, meters
lrm = [-7783653.6376670003, 2223901.0393329998] #lower right, meters

sinu = pyproj.Proj(proj='sinu', lon_0=0.0, x_0=0.0, y_0=0.0)
m = Basemap(projection='sinu', resolution=None, lon_0=0.0)

print "ULM: " + str(ulm)
print "Proj: " + str(sinu(ulm[0], ulm[1], inverse=True))
print "Basemap: " + str(m(ulm[0], ulm[1], inverse=True))
#----- end code --------

This gives:
ULM: [-8895604.1573329996, 3335851.5589999999]
Proj: (-92.327237416031437, 30.141972433747089)
Basemap: (-159.99950210056144, -59.99995206181125)

I'm sure I'm missing something really simple, but I've read a lot of
documentation and I'm not sure what.

Many thanks for any help.

Best,
Tim

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : Jeffrey S. Whitaker: NOAA Physical Sciences Laboratory

_Timothy_W_Hilton · May 4, 2010, 10:25pm

Hi Jeff,

Thanks very much for your response. As you noted, I do not understand
the Basemap global sinusoidal coordinate system. Does this statement
not set up a global sinusoidal cartesian coordinate system centered at
(lon = 0.0, lat = 0.0)?

m = Basemap(projection='sinu', resolution=None, lon_0=0.0, lat_0=0.0)

If so, I would expect m(0.0, 0.0) to return (0.0, 0.0) and m(0.0, 0.0,
inverse=True) to return (0.0, 0.0). Instead, I get:

m(0.0,0.0)

(20015077.371199999, 10007538.6856)

m(0.0,0.0,inverse=True)

(-176.20919036912957, -89.999999999808395)

Sorry if I am being obtuse. Many thanks for your help.

-Tim

···

On Tue, May 2010, 04 at 04:01:21PM -0600, Jeff Whitaker wrote:

On 5/4/10 2:03 PM, Timothy W. Hilton wrote:
>Hello matplotlib users,
>
>I am having trouble understanding the coordinate transformations in
>Basemap and pyproj. I have gridded MODIS vegetation data, with upper
>left corner and lower right corner given in projection coordinates
>(meters). I want to contour the data with Basemap. The data are in a
>sinusoidal projection, but the coordinates do not correspond to what
>Basemap seems to expect.
>
>The code below illustrates the problem. Proj translates the upper
>left to lat/lon correctly (-92.327237416031437, 30.141972433747089),
>while Basemap does not.
>
>#-------- code --------
>from mpl_toolkits.basemap import Basemap
>from mpl_toolkits.basemap import pyproj
>
>ulm = [-8895604.1573329996, 3335851.5589999999] #upper left, meters
>lrm = [-7783653.6376670003, 2223901.0393329998] #lower right, meters
>
>sinu = pyproj.Proj(proj='sinu', lon_0=0.0, x_0=0.0, y_0=0.0)
>m = Basemap(projection='sinu', resolution=None, lon_0=0.0)
>
>print "ULM: " + str(ulm)
>print "Proj: " + str(sinu(ulm[0], ulm[1], inverse=True))
>print "Basemap: " + str(m(ulm[0], ulm[1], inverse=True))
>#----- end code --------
>
>This gives:
>ULM: [-8895604.1573329996, 3335851.5589999999]
>Proj: (-92.327237416031437, 30.141972433747089)
>Basemap: (-159.99950210056144, -59.99995206181125)
>
>I'm sure I'm missing something really simple, but I've read a lot of
>documentation and I'm not sure what.
>
>Many thanks for any help.
>
>Best,
>Tim

Tim: Basemap is using pyproj under the hood, but only supports a
subset of possible proj4 projections. The basemap sinusoidal
projection is global - you can't specify a subregion of the globe.
I think that's where the discrepancy is coming from. I'm sure
there's a way to plot your MODIS data on a global sinusoidal
projection - but it will involve transforming the coordinates to the
Basemap global sinuosidal coordinate system.

-Jeff

_Jeff_Whitaker1 · May 4, 2010, 10:36pm

Hi Jeff,

Thanks very much for your response. As you noted, I do not understand
the Basemap global sinusoidal coordinate system. Does this statement
not set up a global sinusoidal cartesian coordinate system centered at
(lon = 0.0, lat = 0.0)?

Tim: It's a global sinusoidal projection, but x=0,y=0 is not at lon_0,lat_0.

>> from mpl_toolkits.basemap import Basemap
>> m = Basemap(projection='sinu', resolution=None, lon_0=0.0, lat_0=0.0)
>> print m(0,0,inverse=True)
(-176.20919036912957, -89.999999999808395)

I forget why I did it this way, but I think it has to do with the fact that the matplotlib coordinate system has 0,0 in the lower left corner, not the middle.

At any rate, apply a offset to x and y to map to your global coordinate system.

-Jeff

···

On 5/4/10 4:25 PM, Timothy W. Hilton wrote:

m = Basemap(projection='sinu', resolution=None, lon_0=0.0, lat_0=0.0)

If so, I would expect m(0.0, 0.0) to return (0.0, 0.0) and m(0.0, 0.0,
inverse=True) to return (0.0, 0.0). Instead, I get:


m(0.0,0.0)


(20015077.371199999, 10007538.6856)


m(0.0,0.0,inverse=True)


(-176.20919036912957, -89.999999999808395)

Sorry if I am being obtuse. Many thanks for your help.

-Tim

On Tue, May 2010, 04 at 04:01:21PM -0600, Jeff Whitaker wrote:


On 5/4/10 2:03 PM, Timothy W. Hilton wrote:


Hello matplotlib users,

I am having trouble understanding the coordinate transformations in
Basemap and pyproj. I have gridded MODIS vegetation data, with upper
left corner and lower right corner given in projection coordinates
(meters). I want to contour the data with Basemap. The data are in a
sinusoidal projection, but the coordinates do not correspond to what
Basemap seems to expect.

The code below illustrates the problem. Proj translates the upper
left to lat/lon correctly (-92.327237416031437, 30.141972433747089),
while Basemap does not.

#-------- code --------


>from mpl_toolkits.basemap import Basemap
>from mpl_toolkits.basemap import pyproj


ulm = [-8895604.1573329996, 3335851.5589999999] #upper left, meters
lrm = [-7783653.6376670003, 2223901.0393329998] #lower right, meters

sinu = pyproj.Proj(proj='sinu', lon_0=0.0, x_0=0.0, y_0=0.0)
m = Basemap(projection='sinu', resolution=None, lon_0=0.0)

print "ULM: " + str(ulm)
print "Proj: " + str(sinu(ulm[0], ulm[1], inverse=True))
print "Basemap: " + str(m(ulm[0], ulm[1], inverse=True))
#----- end code --------

This gives:
ULM: [-8895604.1573329996, 3335851.5589999999]
Proj: (-92.327237416031437, 30.141972433747089)
Basemap: (-159.99950210056144, -59.99995206181125)

I'm sure I'm missing something really simple, but I've read a lot of
documentation and I'm not sure what.

Many thanks for any help.

Best,
Tim


Tim: Basemap is using pyproj under the hood, but only supports a
subset of possible proj4 projections. The basemap sinusoidal
projection is global - you can't specify a subregion of the globe.
I think that's where the discrepancy is coming from. I'm sure
there's a way to plot your MODIS data on a global sinusoidal
projection - but it will involve transforming the coordinates to the
Basemap global sinuosidal coordinate system.

-Jeff

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : Jeffrey S. Whitaker: NOAA Physical Sciences Laboratory