basemap and omerc

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 <jswhit@...146... > <mailto:jswhit@…146…>> 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 <jswhit@...146... > <mailto:jswhit@…146…> > > <mailto:jswhit@…146…>> 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 <jswhit@...146... > <mailto:jswhit@…146…> > > <mailto:jswhit@…146…> > > > <mailto:jswhit@…146… > <mailto:jswhit@…146…>>> 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 :
    Jeffrey.S.Whitaker@...259... <mailto:Jeffrey.S.Whitaker@…259…>
    > <mailto:Jeffrey.S.Whitaker@…259…
    <mailto:Jeffrey.S.Whitaker@…259…>>
    > 325 Broadway Office : Skaggs Research Cntr 1D-124
    > Boulder, CO, USA 80303-3328 Web : Jeffrey S. Whitaker: NOAA Physical Sciences Laboratory
    >
    ------------------------------------------------------------------------
    >

    --
    Jeffrey S. Whitaker Phone : (303)497-6313
    Meteorologist FAX : (303)497-6449
    NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
    <mailto:Jeffrey.S.Whitaker@…259…>
    325 Broadway Office : Skaggs Research Cntr 1D-124
    Boulder, CO, USA 80303-3328 Web : Jeffrey S. Whitaker: NOAA Physical Sciences Laboratory

--
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-124
Boulder, CO, USA 80303-3328 Web : Jeffrey S. Whitaker: NOAA Physical Sciences Laboratory

That works fine now, thanks very much for your help.

-Evan

···

On Wed, Feb 13, 2008 at 2:16 PM, Jeff Whitaker <jswhit@…146…> 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 <jswhit@…146… > > mailto:jswhit@...146...> 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](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 <jswhit@...146... >  > >     <mailto:jswhit@...146...> > >     > <mailto:jswhit@...146... <mailto:jswhit@...146...>>> 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 <jswhit@...146...
<mailto:jswhit@...146...>
>     <mailto:jswhit@...146... <mailto:jswhit@...146...>>
>     > <mailto:jswhit@...146... <mailto:jswhit@...146...>
<mailto:jswhit@...146... <mailto:jswhit@...146...>>>> 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  :
Jeffrey.S.Whitaker@...259... <mailto:Jeffrey.S.Whitaker@...259...>
>     <mailto:Jeffrey.S.Whitaker@...259...
<mailto:Jeffrey.S.Whitaker@...259...>>
>     325 Broadway                Office : Skaggs Research Cntr 1D-124
>     Boulder, CO, USA 80303-3328 Web    : [http://tinyurl.com/5telg](http://tinyurl.com/5telg)
>
>
>
>
------------------------------------------------------------------------
>
--
Jeffrey S. Whitaker         Phone  : (303)497-6313
Meteorologist               FAX    : (303)497-6449
NOAA/OAR/PSD  R/PSD1        Email  : Jeffrey.S.Whitaker@...259...
<mailto:Jeffrey.S.Whitaker@...259...>
325 Broadway                Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web    : [http://tinyurl.com/5telg](http://tinyurl.com/5telg)

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-124

Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg