From: Anja Roesel <anja.roesel@zm...>  20090806 16:18:38

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 D20146 Hamburg Germany Tel. +49 40 42838 5430 Fax: +49 40 42838 7471 
From: Jeff Whitaker <jswhit@fa...>  20090806 16:26:53

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)4976313 Meteorologist FAX : (303)4976449 NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@... 325 Broadway Office : Skaggs Research Cntr 1D113 Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg 
From: Anja Roesel <anja.roesel@zm...>  20090807 08:54:21

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 > 
From: Anja Roesel <anja.roesel@zm...>  20090807 11:09:35

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 D20146 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 > 
From: Jeff Whitaker <jswhit@fa...>  20090807 12:21:49

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.], > [ 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 