basemap and omerc

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

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?

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.’)

-Evan

···

On Feb 13, 2008 5:14 AM, Jeff Whitaker <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

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

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 <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...> 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…

325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg