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 = 20000.
> > > > 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 = 200000.
> > > 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]>
> > > <mailto:[EMAIL PROTECTED] <mailto:[EMAIL PROTECTED]>>
> > > > <mailto:[EMAIL PROTECTED] <mailto:[EMAIL PROTECTED]>
> > <mailto:[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 = 20000.
> > > > > 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/PSD1 FAX : (303)497-6449
> > > > 325 Broadway Boulder, CO, USA 80305-3328
> > > >
> > > >
> > >
> > >
> > > --
> > > Jeffrey S. Whitaker Phone : (303)497-6313
> > > Meteorologist FAX : (303)497-6449
> > > NOAA/OAR/PSD R/PSD1 Email :
> > [EMAIL PROTECTED] <mailto:[EMAIL PROTECTED]>
> > > <mailto:[EMAIL PROTECTED]
> > <mailto:[EMAIL PROTECTED]>>
> > > 325 Broadway Office : Skaggs Research Cntr
> 1D-124
> > > Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
> > >
> > >
> > >
> > >
> >
> ------------------------------------------------------------------------
> > >
> >
> >
> > --
> > Jeffrey S. Whitaker Phone : (303)497-6313
> > Meteorologist FAX : (303)497-6449
> > NOAA/OAR/PSD R/PSD1 Email : [EMAIL PROTECTED]
> > <mailto:[EMAIL PROTECTED]>
> > 325 Broadway Office : Skaggs Research Cntr 1D-124
> > Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
> >
> >
>
>
> --
> Jeffrey S. Whitaker Phone : (303)497-6313
> Meteorologist FAX : (303)497-6449
> NOAA/OAR/PSD R/PSD1 Email : [EMAIL PROTECTED]
> 325 Broadway Office : Skaggs Research Cntr 1D-124
> Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
>
>
That works fine now, thanks very much for your help.
-Evan
-------------------------------------------------------------------------
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse0120000070mrt/direct/01/
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
