Basemap and UTM Coordinates

Anja_Roesel · August 6, 2009, 3:59pm

Dear list,

I am a matplotlib newbie and I have some simple problems with the
coordinate reprojection.

I have a landsat scene in UTM Projection and I would like to plot it in
a polarstereograhic projection (it is in the Arctic)

I tried it like this:
m=Basemap(resolution='i',projection='npstere',lon_0=-45,boundinglat=70)
XP,YP = m(X_utm,Y_utm)

but the outcome is nonsense (an array with 1.0's)
I am wondering if there is an option to specify the input parameters
like utm zone etc.

Thanks for any help and ideas!
Anja

···

--
Anja Rösel
Center for Marine and Atmospheric Research
Institute of Oceanography, University of Hamburg
Bundesstrasse 53
D-20146 Hamburg
Germany
Tel. +49 40 42838 5430 Fax: +49 40 42838 7471

_Jeff_Whitaker1 · August 6, 2009, 4:26pm

Anja Roesel wrote:

Dear list,

I am a matplotlib newbie and I have some simple problems with the
coordinate reprojection.

I have a landsat scene in UTM Projection and I would like to plot it in
a polarstereograhic projection (it is in the Arctic)

I tried it like this:
m=Basemap(resolution='i',projection='npstere',lon_0=-45,boundinglat=70)
XP,YP = m(X_utm,Y_utm)

but the outcome is nonsense (an array with 1.0's)
I am wondering if there is an option to specify the input parameters
like utm zone etc.

Thanks for any help and ideas!
Anja

Anja: You need to pass the latitude and longitude values (in degrees) to the Basemap instance when converting to projection coordinates. So, to convert from UTM coordinates to polar stereographic coordinates you will need to do something like this:

map1 = Basemap(<parameters for transverse mercator projection>)

lons, lats = map1(x, y, inverse=True) # x and y are projection coordinates on original UTM grid.

# lons and lats are now latitudes and longitudes of UTM grid (in degrees)

map2 =

Basemap(resolution='i',projection='npstere',lon_0=-45,boundinglat=70)

x, y = map2(lons, lats)

# x,y are now polar stereo coordinates of UTM grid.

Or, if you already have the latitudes and logintudes of the original UTM grid you can skip the first two lines and just pass those to the stereographic Basemap instance.

-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

Anja_Roesel · August 7, 2009, 8:54am

Jeff,

great, thank you!

···

--
Anja Rösel

Jeff Whitaker wrote:

Anja Roesel wrote:

Dear list,

I am a matplotlib newbie and I have some simple problems with the
coordinate reprojection.

I have a landsat scene in UTM Projection and I would like to plot it in
a polarstereograhic projection (it is in the Arctic)

I tried it like this:
m=Basemap(resolution='i',projection='npstere',lon_0=-45,boundinglat=70)
XP,YP = m(X_utm,Y_utm)

but the outcome is nonsense (an array with 1.0's)
I am wondering if there is an option to specify the input parameters
like utm zone etc.

Thanks for any help and ideas!
Anja

Anja: You need to pass the latitude and longitude values (in degrees)
to the Basemap instance when converting to projection coordinates. So,
to convert from UTM coordinates to polar stereographic coordinates you
will need to do something like this:

map1 = Basemap(<parameters for transverse mercator projection>)

lons, lats = map1(x, y, inverse=True) # x and y are projection
coordinates on original UTM grid.

# lons and lats are now latitudes and longitudes of UTM grid (in degrees)

map2 =

Basemap(resolution='i',projection='npstere',lon_0=-45,boundinglat=70)

x, y = map2(lons, lats)

# x,y are now polar stereo coordinates of UTM grid.

Or, if you already have the latitudes and logintudes of the original UTM
grid you can skip the first two lines and just pass those to the
stereographic Basemap instance.

-Jeff

Anja_Roesel · August 7, 2009, 11:09am

Thanks Jeff, that's actually what I wanted to know.
Unfortunately I guess, there is a small difference between Transverse
Mercator and Universal Transverse Mercator Projections
Do I have the possibility to define somewere the UTM zone and the Datum
of my Projection.
I tried it like Jeff suggested,

map1 = Basemap(projection='tmerc',
resolution='i',llcrnrlon=-128.3841365,llcrnrlat=68.3952899,urcrnrlon=-121.5846218,urcrnrlat=70.7793990,
lon_0=-124.98437915,lat_0=69.587344449999989)

lons,lats = map1(X_utm,Y_utm,inverse='True')

but the lons and lats
look like this:

In [53]: lons
Out[53]:
array([[-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       ...,
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915]])

In [54]: lats
Out[54]:
array([[ 90., 90., 90., ..., 90., 90., 90.],
       [ 90., 90., 90., ..., 90., 90., 90.],
       [ 90., 90., 90., ..., 90., 90., 90.],
       ...,
       [ 90., 90., 90., ..., 90., 90., 90.],
       [ 90., 90., 90., ..., 90., 90., 90.],
       [ 90., 90., 90., ..., 90., 90., 90.]])

Has anyone an idea?

···

--
Anja Roesel
Center for Marine and Atmospheric Research
Institute of Oceanography, University of Hamburg
Bundesstrasse 53
D-20146 Hamburg
Germany
Tel. +49 40 42838 5430 Fax: +49 40 42838 7471

Anja: You need to pass the latitude and longitude values (in degrees)
to the Basemap instance when converting to projection coordinates. So,
to convert from UTM coordinates to polar stereographic coordinates you
will need to do something like this:

map1 = Basemap(<parameters for transverse mercator projection>)

lons, lats = map1(x, y, inverse=True) # x and y are projection
coordinates on original UTM grid.

# lons and lats are now latitudes and longitudes of UTM grid (in degrees)

map2 =

Basemap(resolution='i',projection='npstere',lon_0=-45,boundinglat=70)

x, y = map2(lons, lats)

# x,y are now polar stereo coordinates of UTM grid.

Or, if you already have the latitudes and logintudes of the original UTM
grid you can skip the first two lines and just pass those to the
stereographic Basemap instance.

-Jeff

_Jeff_Whitaker1 · August 7, 2009, 12:20pm

Anja Roesel wrote:

Thanks Jeff, that's actually what I wanted to know.
Unfortunately I guess, there is a small difference between Transverse
Mercator and Universal Transverse Mercator Projections
Do I have the possibility to define somewere the UTM zone and the Datum
of my Projection.
I tried it like Jeff suggested,

map1 = Basemap(projection='tmerc',
resolution='i',llcrnrlon=-128.3841365,llcrnrlat=68.3952899,urcrnrlon=-121.5846218,urcrnrlat=70.7793990,
lon_0=-124.98437915,lat_0=69.587344449999989)

lons,lats = map1(X_utm,Y_utm,inverse='True')

but the lons and lats
look like this:

In [53]: lons
Out[53]:
array([[-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       ...,
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915],
       [-124.98437915, -124.98437915, -124.98437915, ..., -124.98437915,
        -124.98437915, -124.98437915]])

In [54]: lats
Out[54]:
array([[ 90., 90., 90., ..., 90., 90., 90.],
       [ 90., 90., 90., ..., 90., 90., 90.],
       ...,
       [ 90., 90., 90., ..., 90., 90., 90.],
       [ 90., 90., 90., ..., 90., 90., 90.]])

Has anyone an idea?

Anja: Basemap measures X and Y in meters from the lower left corner of the plot. I don't know how your X_utm and Y_utm are defined, or what units they are. You also may have to specify the major and minor radii of the earth to exactly match the UTM projection definition you are using. Basemap cannot define a UTM projection, the the pyproj module inside basemap does. If you know the UTM zone you could do, for example:

from mpl_toolkits.basemap import pyproj
p = pyproj.Proj(proj='utm',zone=10,ellps='WGS84')
lons, lats = p(X_utm,Y_utm,inverse=True)

you would still need to be careful about the units of X_utm and Y_utm and the ellipse definition, however.

-Jeff