Re: [Matplotlib-users] basemap and omerc

2008-02-14 Thread Evan Mason
On Wed, Feb 13, 2008 at 2:16 PM, Jeff Whitaker [EMAIL PROTECTED] wrote:

 Evan Mason wrote:
  Hi Jeff
 
  Here are the corners:
 
  lon_corners = N.array([-4.09300764,-35.76003475,-43.72330207,
  -12.05627497])
  lat_corners = N.array([41.90278813, 49.2136974, 14.7209971, 7.41008784])
 
  The reason for the differences is that the matlab script is very
  fiddly, lots of trial and error to get the grid in the right place.
  The attraction of using basemap is it allows me to specify the
  corners, so that I have it right first time, that's the idea anyway.
 
  That would be great if you could turn off that rotation, maybe a
  keyword True/False
 
  Thanks, Evan

 Evan:  I've changed Basemap in svn so you can set 'no_rot=True' when
 creating a Basemap instance for the 'omerc' projection to get what you
 want.  If you don't feel like upgrading (since that requires upgrading
 matplotlib to svn head at the same time), this will work in the version
 you have:

 from matplotlib.toolkits.basemap import Basemap, pyproj
 from pylab import *
 p = pyproj.Proj(lon_2=-27.8,lon_1=-19.9,no_rot=True,proj='omerc',\
lat_2=11.0,lat_1=45.5)
 xmax,ymax = p(-4.093,41.9027)  # upper right corner
 xmin,ymin = p(-43.723,14.721)  # lower left corner
 x = linspace(xmin,xmax,35)
 y = linspace(ymin,ymax,35)
 x, y = meshgrid(x,y)
 lonr,latr = p(x,y, inverse=True)
 m = Basemap(llcrnrlon=-60,llcrnrlat=5,\
urcrnrlon=15,urcrnrlat=60,resolution='i')
 m.drawcoastlines()
 m.scatter(lonr.flatten(),latr.flatten(),5,marker='o')
 m.drawmeridians(arange(-60,21,10),labels=[0,0,0,1])
 m.drawparallels(arange(0,61,10),labels=[1,0,0,0])
 show()

 Let me know if this fixes it for you.

 -Jeff
 
 
 
  On Feb 13, 2008 12:56 PM, Jeff Whitaker [EMAIL PROTECTED]
  mailto:[EMAIL PROTECTED] wrote:
 
  Evan Mason wrote:
   Hi Jeff
  
   By losing the memory I mean that the grid is no longer rotated;
 that
   the rotation I introduced through lat1, lon1, lat2, lon2 is
  lost.  If
   you look at the latitude of the two bottom corners you see that
 they
   are the same, they should be different.  In other words, I want my
   great circle not to be the equator or a meridian, instead I want
  it to
   be between lat1, lon1, lat2, lon2.  See for example:
  
 
 http://erg.usgs.gov/isb/pubs/MapProjections/projections.html#mercator
  
   Attached is a png from the matlab script.  Here you can see the
   rotation that I am looking for.  The latitude of the two bottom
   corners is different, unlike what happens presently with my
 basemap
   script.
  
   Thanks, Evan
 
  Evan:  OK, I was confused by your use of the term 'losing the
 memory'.
  Basemap didn't lose the rotation, it never had it in the first
 place.
  It looks like matlab and Basemap define the projection regions
  differently.  They both are right, but are showing you different
  regions
  of the same projection.  The difference is that proj4 (and therefore
  Basemap) automatically rotates the y axis to lie along true north.
  I
  think I know how to modify Basemap to display the region you want,
 by
  turning off that rotation.  Can you send me the lat/lon values of
  the 4
  corners of the region that matlab produces?
 
  -Jeff
 
  P.S. I don't know if this is relevant or not, but you appear to be
  giving matlab different points to define the center of the
 projection
  than you did in Basemap (the lons you gave matlab are
  -23.75,-28.25, the
  lons you give in Basemap are -27.8 and 19.9).
  
  
  
   On Feb 13, 2008 10:48 AM, Jeff Whitaker [EMAIL PROTECTED]
  mailto:[EMAIL PROTECTED]
   mailto:[EMAIL PROTECTED] mailto:[EMAIL PROTECTED] wrote:
  
   Evan Mason wrote:
Thanks for the replies.  The map you produced, Jeff, looks
  as it
should.  However, I am trying to make an ocean model grid,
  and so I
require two 2d arrays of lon and lat, at my desired grid
  spacing.
This is why I try the steps:
   
dl = 2.
nx = int((M.xmax - M.xmin) / dl) + 1
ny = int((M.ymax - M.ymin) / dl) + 1
lonr, latr = M.makegrid(nx, ny)   - it seems to be here
  that it
   loses
'memory' of omerc projection that I specified, and maybe
  there is a
bug here?
  
   Evan:  Why do you say it 'loses' memory of the projection?
   The values
   look fine to me - they are just equally spaced points in map
   projection
   coordinates that cover the map projection region.  Take a
  look at
  
   M = Basemap(projection = 'omerc',   \
   resolution  = 'l',   \
  llcrnrlon  = -43.7,   \
  llcrnrlat   = 14.7,\
 

Re: [Matplotlib-users] basemap and omerc

2008-02-13 Thread Jeff Whitaker
Evan Mason wrote:
 Hi, I am having some problems using the oblique mercator projection in 
 basemap.  I want to define a rectangular orthogonal grid, rotated 
 clockwise by about 13 degrees.  I want to define grid cells of size, 
 say, about 20x20 km.  The script I have so far is below.  The problem 
 is that at some point (the makegrid step) I lose the rotation, as seen 
 in the plot.

 I'd appreciate any help with this, thanks, Evan


 from matplotlib.toolkits.basemap import Basemap

 M = Basemap(projection = 'omerc',   \
resolution  = None,   \
llcrnrlon  = -43.7,   \
llcrnrlat   = 14.7,\
urcrnrlon = -4.0,\
urcrnrlat  = 41.9,\
lat_2   = 11.0,\
lat_1   = 45.5,\
lon_2  = -27.8,   \
lon_1  = -19.9)

 dl = 2.
 nx = int((M.xmax - M.xmin) / dl) + 1
 ny = int((M.ymax - M.ymin) / dl) + 1

 lonr, latr = M.makegrid(nx, ny)

 plot(lonr, latr, 'c.')
 show()

Evan:  I have to admit, I'm not too familiar with the Oblique Mercator 
projection.  What exactly should it look like?

If I plot

M = Basemap(projection = 'omerc',   \
   resolution  = 'l',   \
   llcrnrlon  = -43.7,   \
   llcrnrlat   = 14.7,\
   urcrnrlon = -4.0,\
   urcrnrlat  = 41.9,\
   lat_2   = 11.0,\
   lat_1   = 45.5,\
   lon_2  = -27.8,   \
   lon_1  = -19.9)
M.drawcoastlines()
M.drawparallels(arange(10,51,10))
M.drawmeridians(arange(-50,1,10))
M.show()

I see a reasonable looking map, but then I don't really know exactly 
what to expect.

It seems that there are two ways to specify oblique mercator in proj4

1) by specifying 2 points (lon1,lat1), (lon2,lat2) along the central line
2) by specifying a central point and an azimuth that passes through the 
central point.

Basemap uses (1), but every example on the web I've seen uses (2).  It 
could be there are bugs in (1), and (2) would produce more reasonable 
results in your case.  If you can give me an example of what your map 
*should* look like,  it would help a lot.

-Jeff




-- 
Jeffrey S. Whitaker Phone : (303)497-6313
NOAA/OAR/CDC  R/PSD1FAX   : (303)497-6449
325 BroadwayBoulder, CO, USA 80305-3328


-
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse012070mrt/direct/01/
___
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users


[Matplotlib-users] basemap and omerc

2008-02-13 Thread Tommy Grav

On Feb 12, 2008, at 8:39 PM, Evan Mason wrote:

 dl = 2.
 nx = int((M.xmax - M.xmin) / dl) + 1
 ny = int((M.ymax - M.ymin) / dl) + 1

 lonr, latr = M.makegrid(nx, ny)

 plot(lonr, latr, 'c.')
 show()

I think you might be looking for M.plot() rather than plot()??
plot() will just overwrite the basemap in the figure (although
that might be possible to fix with some axis-magic, but that is
far beyond my knowledge of matplotlib).

Cheers
   Tommy


-
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse012070mrt/direct/01/
___
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users


Re: [Matplotlib-users] basemap and omerc

2008-02-13 Thread Jeff Whitaker
Evan Mason wrote:
 Thanks for the replies.  The map you produced, Jeff, looks as it 
 should.  However, I am trying to make an ocean model grid, and so I 
 require two 2d arrays of lon and lat, at my desired grid spacing.  
 This is why I try the steps:

 dl = 2.
 nx = int((M.xmax - M.xmin) / dl) + 1
 ny = int((M.ymax - M.ymin) / dl) + 1
 lonr, latr = M.makegrid(nx, ny)   - it seems to be here that it loses 
 'memory' of omerc projection that I specified, and maybe there is a 
 bug here?

Evan:  Why do you say it 'loses' memory of the projection?  The values 
look fine to me - they are just equally spaced points in map projection 
coordinates that cover the map projection region.  Take a look at

M = Basemap(projection = 'omerc',   \
 resolution  = 'l',   \
llcrnrlon  = -43.7,   \
llcrnrlat   = 14.7,\
urcrnrlon = -4.0,\
urcrnrlat  = 41.9,\
lat_2   = 11.0,\
lat_1   = 45.5,\
lon_2  = -27.8,   \
lon_1  = -19.9)
dl = 20.
nx = int((M.xmax - M.xmin) / dl) + 1
ny = int((M.ymax - M.ymin) / dl) + 1
lonr, latr,x,y= M.makegrid(nx, ny, returnxy=True)
M.drawcoastlines()
M.scatter(x.flatten(), y.flatten(),5,marker='o')
M.drawparallels(arange(10,51,10))
M.drawmeridians(arange(-50,1,10))
show()

 If you have matlab, the following lines do what I am looking for:

 incx = 0.00310/2;
 incy = 0.00306/2;
 Xstr = -0.275;
 Xend = 0.275;
 Ystr  = 0.17;
 Yend = 0.8;
 X = [Xstr:incx:Xend];
 Y = [Ystr:incy:Yend];
 [XX,YY]= meshgrid(X,Y);
 [Lonr,Latr] = m_xy2ll(XX,YY);
 m_proj('Oblique Mercator','lon',[-23.75 -28.25],'lat',[45.5 
 11],'direction','vertical');
 plot(Lonr, Latr, 'c.')

Sorry, I don't have matlab - but it looks at first glance like it ought 
to be doing the same thing.

-Jeff



 -Evan





 On Feb 13, 2008 5:14 AM, Jeff Whitaker [EMAIL PROTECTED] 
 mailto:[EMAIL PROTECTED] wrote:

 Evan Mason wrote:
  Hi, I am having some problems using the oblique mercator
 projection in
  basemap.  I want to define a rectangular orthogonal grid, rotated
  clockwise by about 13 degrees.  I want to define grid cells of size,
  say, about 20x20 km.  The script I have so far is below.  The
 problem
  is that at some point (the makegrid step) I lose the rotation,
 as seen
  in the plot.
 
  I'd appreciate any help with this, thanks, Evan
 
 
  from matplotlib.toolkits.basemap import Basemap
 
  M = Basemap(projection = 'omerc',   \
 resolution  = None,   \
 llcrnrlon  = -43.7,   \
 llcrnrlat   = 14.7,\
 urcrnrlon = -4.0,\
 urcrnrlat  = 41.9,\
 lat_2   = 11.0,\
 lat_1   = 45.5,\
 lon_2  = -27.8,   \
 lon_1  = -19.9)
 
  dl = 2.
  nx = int((M.xmax - M.xmin) / dl) + 1
  ny = int((M.ymax - M.ymin) / dl) + 1
 
  lonr, latr = M.makegrid(nx, ny)
 
  plot(lonr, latr, 'c.')
  show()

 Evan:  I have to admit, I'm not too familiar with the Oblique Mercator
 projection.  What exactly should it look like?

 If I plot

 M = Basemap(projection = 'omerc',   \
   resolution  = 'l',   \
   llcrnrlon  = -43.7,   \
   llcrnrlat   = 14.7,\
   urcrnrlon = -4.0,\
   urcrnrlat  = 41.9,\
   lat_2   = 11.0,\
   lat_1   = 45.5,\
   lon_2  = -27.8,   \
   lon_1  = -19.9)
 M.drawcoastlines()
 M.drawparallels(arange(10,51,10))
 M.drawmeridians(arange(-50,1,10))
 M.show()

 I see a reasonable looking map, but then I don't really know exactly
 what to expect.

 It seems that there are two ways to specify oblique mercator in proj4

 1) by specifying 2 points (lon1,lat1), (lon2,lat2) along the
 central line
 2) by specifying a central point and an azimuth that passes
 through the
 central point.

 Basemap uses (1), but every example on the web I've seen uses (2).  It
 could be there are bugs in (1), and (2) would produce more reasonable
 results in your case.  If you can give me an example of what your map
 *should* look like,  it would help a lot.

 -Jeff




 --
 Jeffrey S. Whitaker Phone : (303)497-6313
 NOAA/OAR/CDC  R/PSD1FAX   : (303)497-6449
 325 BroadwayBoulder, CO, USA 80305-3328




-- 
Jeffrey S. Whitaker Phone  : (303)497-6313
Meteorologist   FAX: (303)497-6449
NOAA/OAR/PSD  R/PSD1Email  : [EMAIL PROTECTED]
325 BroadwayOffice : Skaggs Research Cntr 1D-124

Re: [Matplotlib-users] basemap and omerc

2008-02-13 Thread Evan Mason
Thanks for the replies.  The map you produced, Jeff, looks as it should.
However, I am trying to make an ocean model grid, and so I require two 2d
arrays of lon and lat, at my desired grid spacing.  This is why I try the
steps:

dl = 2.
nx = int((M.xmax - M.xmin) / dl) + 1
ny = int((M.ymax - M.ymin) / dl) + 1
lonr, latr = M.makegrid(nx, ny)   - it seems to be here that it loses
'memory' of omerc projection that I specified, and maybe there is a bug
here?

If you have matlab, the following lines do what I am looking for:

incx = 0.00310/2;
incy = 0.00306/2;
Xstr = -0.275;
Xend = 0.275;
Ystr  = 0.17;
Yend = 0.8;
X = [Xstr:incx:Xend];
Y = [Ystr:incy:Yend];
[XX,YY]= meshgrid(X,Y);
[Lonr,Latr] = m_xy2ll(XX,YY);
m_proj('Oblique Mercator','lon',[-23.75
-28.25],'lat',[45.511],'direction','vertical');
plot(Lonr, Latr, 'c.')



-Evan





On Feb 13, 2008 5:14 AM, Jeff Whitaker [EMAIL PROTECTED] wrote:

 Evan Mason wrote:
  Hi, I am having some problems using the oblique mercator projection in
  basemap.  I want to define a rectangular orthogonal grid, rotated
  clockwise by about 13 degrees.  I want to define grid cells of size,
  say, about 20x20 km.  The script I have so far is below.  The problem
  is that at some point (the makegrid step) I lose the rotation, as seen
  in the plot.
 
  I'd appreciate any help with this, thanks, Evan
 
 
  from matplotlib.toolkits.basemap import Basemap
 
  M = Basemap(projection = 'omerc',   \
 resolution  = None,   \
 llcrnrlon  = -43.7,   \
 llcrnrlat   = 14.7,\
 urcrnrlon = -4.0,\
 urcrnrlat  = 41.9,\
 lat_2   = 11.0,\
 lat_1   = 45.5,\
 lon_2  = -27.8,   \
 lon_1  = -19.9)
 
  dl = 2.
  nx = int((M.xmax - M.xmin) / dl) + 1
  ny = int((M.ymax - M.ymin) / dl) + 1
 
  lonr, latr = M.makegrid(nx, ny)
 
  plot(lonr, latr, 'c.')
  show()

 Evan:  I have to admit, I'm not too familiar with the Oblique Mercator
 projection.  What exactly should it look like?

 If I plot

 M = Basemap(projection = 'omerc',   \
   resolution  = 'l',   \
   llcrnrlon  = -43.7,   \
   llcrnrlat   = 14.7,\
   urcrnrlon = -4.0,\
   urcrnrlat  = 41.9,\
   lat_2   = 11.0,\
   lat_1   = 45.5,\
   lon_2  = -27.8,   \
   lon_1  = -19.9)
 M.drawcoastlines()
 M.drawparallels(arange(10,51,10))
 M.drawmeridians(arange(-50,1,10))
 M.show()

 I see a reasonable looking map, but then I don't really know exactly
 what to expect.

 It seems that there are two ways to specify oblique mercator in proj4

 1) by specifying 2 points (lon1,lat1), (lon2,lat2) along the central line
 2) by specifying a central point and an azimuth that passes through the
 central point.

 Basemap uses (1), but every example on the web I've seen uses (2).  It
 could be there are bugs in (1), and (2) would produce more reasonable
 results in your case.  If you can give me an example of what your map
 *should* look like,  it would help a lot.

 -Jeff




 --
 Jeffrey S. Whitaker Phone : (303)497-6313
 NOAA/OAR/CDC  R/PSD1FAX   : (303)497-6449
 325 BroadwayBoulder, CO, USA 80305-3328


-
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse012070mrt/direct/01/___
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users


Re: [Matplotlib-users] basemap and omerc

2008-02-13 Thread Evan Mason
By losing the memory I mean that the grid is no longer rotated; that the
rotation I introduced through lat1, lon1, lat2, lon2 is lost.  If you look
at the latitude of the two bottom corners you see that they are the same,
they should be different - for the matlab script they are different.  In
other words, I want my great circle not to be the equator or a meridian,
instead I want it to be between lat1, lon1, lat2, lon2.  See for example:
http://erg.usgs.gov/isb/pubs/MapProjections/projections.html#mercator

At present, basemap seems to be reverting to a standard mercator projection.

-Evan




On Feb 13, 2008 10:48 AM, Jeff Whitaker [EMAIL PROTECTED] wrote:

 Evan Mason wrote:
  Thanks for the replies.  The map you produced, Jeff, looks as it
  should.  However, I am trying to make an ocean model grid, and so I
  require two 2d arrays of lon and lat, at my desired grid spacing.
  This is why I try the steps:
 
  dl = 2.
  nx = int((M.xmax - M.xmin) / dl) + 1
  ny = int((M.ymax - M.ymin) / dl) + 1
  lonr, latr = M.makegrid(nx, ny)   - it seems to be here that it loses
  'memory' of omerc projection that I specified, and maybe there is a
  bug here?

 Evan:  Why do you say it 'loses' memory of the projection?  The values
 look fine to me - they are just equally spaced points in map projection
 coordinates that cover the map projection region.  Take a look at

 M = Basemap(projection = 'omerc',   \
 resolution  = 'l',   \
llcrnrlon  = -43.7,   \
llcrnrlat   = 14.7,\
urcrnrlon = -4.0,\
urcrnrlat  = 41.9,\
lat_2   = 11.0,\
lat_1   = 45.5,\
lon_2  = -27.8,   \
lon_1  = -19.9)
 dl = 20.
 nx = int((M.xmax - M.xmin) / dl) + 1
 ny = int((M.ymax - M.ymin) / dl) + 1
 lonr, latr,x,y= M.makegrid(nx, ny, returnxy=True)
 M.drawcoastlines()
 M.scatter(x.flatten(), y.flatten(),5,marker='o')
 M.drawparallels(arange(10,51,10))
 M.drawmeridians(arange(-50,1,10))
 show()
 
  If you have matlab, the following lines do what I am looking for:
 
  incx = 0.00310/2;
  incy = 0.00306/2;
  Xstr = -0.275;
  Xend = 0.275;
  Ystr  = 0.17;
  Yend = 0.8;
  X = [Xstr:incx:Xend];
  Y = [Ystr:incy:Yend];
  [XX,YY]= meshgrid(X,Y);
  [Lonr,Latr] = m_xy2ll(XX,YY);
  m_proj('Oblique Mercator','lon',[-23.75 -28.25],'lat',[45.5
  11],'direction','vertical');
  plot(Lonr, Latr, 'c.')

 Sorry, I don't have matlab - but it looks at first glance like it ought
 to be doing the same thing.

 -Jeff
 
 
 
  -Evan
 
 
 
 
 
  On Feb 13, 2008 5:14 AM, Jeff Whitaker [EMAIL PROTECTED]
  mailto:[EMAIL PROTECTED] wrote:
 
  Evan Mason wrote:
   Hi, I am having some problems using the oblique mercator
  projection in
   basemap.  I want to define a rectangular orthogonal grid, rotated
   clockwise by about 13 degrees.  I want to define grid cells of
 size,
   say, about 20x20 km.  The script I have so far is below.  The
  problem
   is that at some point (the makegrid step) I lose the rotation,
  as seen
   in the plot.
  
   I'd appreciate any help with this, thanks, Evan
  
  
   from matplotlib.toolkits.basemap import Basemap
  
   M = Basemap(projection = 'omerc',   \
  resolution  = None,   \
  llcrnrlon  = -43.7,   \
  llcrnrlat   = 14.7,\
  urcrnrlon = -4.0,\
  urcrnrlat  = 41.9,\
  lat_2   = 11.0,\
  lat_1   = 45.5,\
  lon_2  = -27.8,   \
  lon_1  = -19.9)
  
   dl = 2.
   nx = int((M.xmax - M.xmin) / dl) + 1
   ny = int((M.ymax - M.ymin) / dl) + 1
  
   lonr, latr = M.makegrid(nx, ny)
  
   plot(lonr, latr, 'c.')
   show()
 
  Evan:  I have to admit, I'm not too familiar with the Oblique
 Mercator
  projection.  What exactly should it look like?
 
  If I plot
 
  M = Basemap(projection = 'omerc',   \
resolution  = 'l',   \
llcrnrlon  = -43.7,   \
llcrnrlat   = 14.7,\
urcrnrlon = -4.0,\
urcrnrlat  = 41.9,\
lat_2   = 11.0,\
lat_1   = 45.5,\
lon_2  = -27.8,   \
lon_1  = -19.9)
  M.drawcoastlines()
  M.drawparallels(arange(10,51,10))
  M.drawmeridians(arange(-50,1,10))
  M.show()
 
  I see a reasonable looking map, but then I don't really know exactly
  what to expect.
 
  It seems that there are two ways to specify oblique mercator in
 proj4
 
  1) by specifying 2 points (lon1,lat1), (lon2,lat2) along the
  central line
  2) by specifying a 

Re: [Matplotlib-users] basemap and omerc

2008-02-13 Thread Jeff Whitaker
Evan Mason wrote:
 Hi Jeff

 Here are the corners:

 lon_corners = N.array([-4.09300764,-35.76003475,-43.72330207, 
 -12.05627497])
 lat_corners = N.array([41.90278813, 49.2136974, 14.7209971, 7.41008784])

 The reason for the differences is that the matlab script is very 
 fiddly, lots of trial and error to get the grid in the right place.  
 The attraction of using basemap is it allows me to specify the 
 corners, so that I have it right first time, that's the idea anyway.

 That would be great if you could turn off that rotation, maybe a 
 keyword True/False

 Thanks, Evan

Evan:  I've changed Basemap in svn so you can set 'no_rot=True' when 
creating a Basemap instance for the 'omerc' projection to get what you 
want.  If you don't feel like upgrading (since that requires upgrading 
matplotlib to svn head at the same time), this will work in the version 
you have:

from matplotlib.toolkits.basemap import Basemap, pyproj
from pylab import *
p = pyproj.Proj(lon_2=-27.8,lon_1=-19.9,no_rot=True,proj='omerc',\
lat_2=11.0,lat_1=45.5)
xmax,ymax = p(-4.093,41.9027)  # upper right corner
xmin,ymin = p(-43.723,14.721)  # lower left corner
x = linspace(xmin,xmax,35)
y = linspace(ymin,ymax,35)
x, y = meshgrid(x,y)
lonr,latr = p(x,y, inverse=True)
m = Basemap(llcrnrlon=-60,llcrnrlat=5,\
urcrnrlon=15,urcrnrlat=60,resolution='i')
m.drawcoastlines()
m.scatter(lonr.flatten(),latr.flatten(),5,marker='o')
m.drawmeridians(arange(-60,21,10),labels=[0,0,0,1])
m.drawparallels(arange(0,61,10),labels=[1,0,0,0])
show()

Let me know if this fixes it for you.

-Jeff



 On Feb 13, 2008 12:56 PM, Jeff Whitaker [EMAIL PROTECTED] 
 mailto:[EMAIL PROTECTED] wrote:

 Evan Mason wrote:
  Hi Jeff
 
  By losing the memory I mean that the grid is no longer rotated; that
  the rotation I introduced through lat1, lon1, lat2, lon2 is
 lost.  If
  you look at the latitude of the two bottom corners you see that they
  are the same, they should be different.  In other words, I want my
  great circle not to be the equator or a meridian, instead I want
 it to
  be between lat1, lon1, lat2, lon2.  See for example:
 
 http://erg.usgs.gov/isb/pubs/MapProjections/projections.html#mercator
 
  Attached is a png from the matlab script.  Here you can see the
  rotation that I am looking for.  The latitude of the two bottom
  corners is different, unlike what happens presently with my basemap
  script.
 
  Thanks, Evan

 Evan:  OK, I was confused by your use of the term 'losing the memory'.
 Basemap didn't lose the rotation, it never had it in the first place.
 It looks like matlab and Basemap define the projection regions
 differently.  They both are right, but are showing you different
 regions
 of the same projection.  The difference is that proj4 (and therefore
 Basemap) automatically rotates the y axis to lie along true north.  I
 think I know how to modify Basemap to display the region you want, by
 turning off that rotation.  Can you send me the lat/lon values of
 the 4
 corners of the region that matlab produces?

 -Jeff

 P.S. I don't know if this is relevant or not, but you appear to be
 giving matlab different points to define the center of the projection
 than you did in Basemap (the lons you gave matlab are
 -23.75,-28.25, the
 lons you give in Basemap are -27.8 and 19.9).
 
 
 
  On Feb 13, 2008 10:48 AM, Jeff Whitaker [EMAIL PROTECTED]
 mailto:[EMAIL PROTECTED]
  mailto:[EMAIL PROTECTED] mailto:[EMAIL PROTECTED] wrote:
 
  Evan Mason wrote:
   Thanks for the replies.  The map you produced, Jeff, looks
 as it
   should.  However, I am trying to make an ocean model grid,
 and so I
   require two 2d arrays of lon and lat, at my desired grid
 spacing.
   This is why I try the steps:
  
   dl = 2.
   nx = int((M.xmax - M.xmin) / dl) + 1
   ny = int((M.ymax - M.ymin) / dl) + 1
   lonr, latr = M.makegrid(nx, ny)   - it seems to be here
 that it
  loses
   'memory' of omerc projection that I specified, and maybe
 there is a
   bug here?
 
  Evan:  Why do you say it 'loses' memory of the projection?
  The values
  look fine to me - they are just equally spaced points in map
  projection
  coordinates that cover the map projection region.  Take a
 look at
 
  M = Basemap(projection = 'omerc',   \
  resolution  = 'l',   \
 llcrnrlon  = -43.7,   \
 llcrnrlat   = 14.7,\
 urcrnrlon = -4.0,\
 urcrnrlat  = 41.9,\
 lat_2   = 11.0,\
 

[Matplotlib-users] basemap and omerc

2008-02-12 Thread Evan Mason
Hi, I am having some problems using the oblique mercator projection in
basemap.  I want to define a rectangular orthogonal grid, rotated clockwise
by about 13 degrees.  I want to define grid cells of size, say, about 20x20
km.  The script I have so far is below.  The problem is that at some point
(the makegrid step) I lose the rotation, as seen in the plot.

I'd appreciate any help with this, thanks, Evan


from matplotlib.toolkits.basemap import Basemap

M = Basemap(projection = 'omerc',   \
   resolution  = None,   \
   llcrnrlon  = -43.7,   \
   llcrnrlat   = 14.7,\
   urcrnrlon = -4.0,\
   urcrnrlat  = 41.9,\
   lat_2   = 11.0,\
   lat_1   = 45.5,\
   lon_2  = -27.8,   \
   lon_1  = -19.9)

dl = 2.
nx = int((M.xmax - M.xmin) / dl) + 1
ny = int((M.ymax - M.ymin) / dl) + 1

lonr, latr = M.makegrid(nx, ny)

plot(lonr, latr, 'c.')
show()
-
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse012070mrt/direct/01/___
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users