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 = 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]>> 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]
> 325 Broadway Office : Skaggs Research Cntr 1D-124
> Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
>
>
-------------------------------------------------------------------------
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
