Screenshot instructions:
Windows
Mac
Red Hat Linux
Ubuntu
Click URL instructions:
Rightclick on ad, choose "Copy Link", then paste here →
(This may not be possible with some types of ads)
From: Evan Mason <evanmason@gm...>  20080213 01:39:34
Attachments:
Message as HTML

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() 
From: Tommy Grav <tgrav@ma...>  20080213 13:58:15

On Feb 12, 2008, at 8:39 PM, Evan Mason wrote: > 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() 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 axismagic, but that is far beyond my knowledge of matplotlib). Cheers Tommy 
From: Jeff Whitaker <jswhit@fa...>  20080213 13:14:28

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)4976313 NOAA/OAR/CDC R/PSD1 FAX : (303)4976449 325 Broadway Boulder, CO, USA 803053328 
From: Evan Mason <evanmason@gm...>  20080213 18:39:23
Attachments:
Message as HTML

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.511],'direction','vertical'); plot(Lonr, Latr, 'c.') Evan On Feb 13, 2008 5:14 AM, Jeff Whitaker <jswhit@...> 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)4976313 > NOAA/OAR/CDC R/PSD1 FAX : (303)4976449 > 325 Broadway Boulder, CO, USA 803053328 > > 
From: Jeff Whitaker <jswhit@fa...>  20080213 18:49:09

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@... > <mailto:jswhit@...>> 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)4976313 > NOAA/OAR/CDC R/PSD1 FAX : (303)4976449 > 325 Broadway Boulder, CO, USA 803053328 > >  Jeffrey S. Whitaker Phone : (303)4976313 Meteorologist FAX : (303)4976449 NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@... 325 Broadway Office : Skaggs Research Cntr 1D124 Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg 
From: Evan Mason <evanmason@gm...>  20080213 20:07:07
Attachments:
Message as HTML

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@...> 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@... > > <mailto:jswhit@...>> 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)4976313 > > NOAA/OAR/CDC R/PSD1 FAX : (303)4976449 > > 325 Broadway Boulder, CO, USA 803053328 > > > > > > >  > Jeffrey S. Whitaker Phone : (303)4976313 > Meteorologist FAX : (303)4976449 > NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@... > 325 Broadway Office : Skaggs Research Cntr 1D124 > Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg > > 
From: Jeff Whitaker <jswhit@fa...>  20080213 22:17:09

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@... > <mailto:jswhit@...>> 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@... > <mailto:jswhit@...> > > <mailto:jswhit@... <mailto:jswhit@...>>> 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@... > <mailto:jswhit@...> > > <mailto:jswhit@... <mailto:jswhit@...>> > > > <mailto:jswhit@... <mailto:jswhit@...> > <mailto:jswhit@... <mailto:jswhit@...>>>> 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)4976313 > > > NOAA/OAR/CDC R/PSD1 FAX : (303)4976449 > > > 325 Broadway Boulder, CO, USA 803053328 > > > > > > > > > > > >  > > Jeffrey S. Whitaker Phone : (303)4976313 > > Meteorologist FAX : (303)4976449 > > NOAA/OAR/PSD R/PSD1 Email : > Jeffrey.S.Whitaker@... <mailto:Jeffrey.S.Whitaker@...> > > <mailto:Jeffrey.S.Whitaker@... > <mailto:Jeffrey.S.Whitaker@...>> > > 325 Broadway Office : Skaggs Research Cntr 1D124 > > Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg > > > > > > > > >  > > > > >  > Jeffrey S. Whitaker Phone : (303)4976313 > Meteorologist FAX : (303)4976449 > NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@... > <mailto:Jeffrey.S.Whitaker@...> > 325 Broadway Office : Skaggs Research Cntr 1D124 > Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg > >  Jeffrey S. Whitaker Phone : (303)4976313 Meteorologist FAX : (303)4976449 NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@... 325 Broadway Office : Skaggs Research Cntr 1D124 Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg 
From: Evan Mason <evanmason@gm...>  20080214 17:28:33
Attachments:
Message as HTML

On Wed, Feb 13, 2008 at 2:16 PM, Jeff Whitaker <jswhit@...> 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@... > > <mailto:jswhit@...>> 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@... > > <mailto:jswhit@...> > > > <mailto:jswhit@... <mailto:jswhit@...>>> 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@... > > <mailto:jswhit@...> > > > <mailto:jswhit@... <mailto:jswhit@...>> > > > > <mailto:jswhit@... <mailto:jswhit@...> > > <mailto:jswhit@... <mailto:jswhit@...>>>> 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)4976313 > > > > NOAA/OAR/CDC R/PSD1 FAX : (303)4976449 > > > > 325 Broadway Boulder, CO, USA 803053328 > > > > > > > > > > > > > > > > >  > > > Jeffrey S. Whitaker Phone : (303)4976313 > > > Meteorologist FAX : (303)4976449 > > > NOAA/OAR/PSD R/PSD1 Email : > > Jeffrey.S.Whitaker@... <mailto:Jeffrey.S.Whitaker@...> > > > <mailto:Jeffrey.S.Whitaker@... > > <mailto:Jeffrey.S.Whitaker@...>> > > > 325 Broadway Office : Skaggs Research Cntr > 1D124 > > > Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg > > > > > > > > > > > > > > >  > > > > > > > > >  > > Jeffrey S. Whitaker Phone : (303)4976313 > > Meteorologist FAX : (303)4976449 > > NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@... > > <mailto:Jeffrey.S.Whitaker@...> > > 325 Broadway Office : Skaggs Research Cntr 1D124 > > Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg > > > > > > >  > Jeffrey S. Whitaker Phone : (303)4976313 > Meteorologist FAX : (303)4976449 > NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@... > 325 Broadway Office : Skaggs Research Cntr 1D124 > Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg > > That works fine now, thanks very much for your help. Evan 