Jeff,
Here's a quick snippet. I've looked at the test.py file provided with the
basemap examples. What I am unclear on are the different ways in which nx
and ny are defined. I would like to have this 'automatically' defined, based
solely on variables from my input object.. say for example a netcdf file
that has len and lon dimensions defined.
Below is my crude stab at it, but I am clearly having some problems. I guess
the point is, maybe it's not possible to have a Basemap instance with
extents beyond the imshow object. Then perhaps I need to make sure that when
I set up the Basemap instance, I pass the H.outlon0 to llcrnrlon for
example. But is that necessary?
Thanks!
#!/usr/bin/env python
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
def plot_imshow_custom(H,transform=True ):
"""
function to automagically plot an mxn array of arbitrary lats/lons
"""
data = H.data
print data.shape
m =
Basemap(projection='npstere',boundinglat=10,lon_0=270,resolution='l')
fig = plt.figure()
ax = fig.gca()
print "Preparing to plot %s with dimensions:" % H.name
print "lon0, numx, dx:"
print H.outlon0, H.numxgrid, H.dxout
print "lat0, numy, dy:"
print H.outlat0, H.numygrid, H.dyout
## set up transformations for the data array
## THIS IS WHERE I NEED SOME HELP:
if m.projection not in ['cyl','merc','mill']:
lats = np.arange( H.outlat0, ( H.outlat0 + ( H.numygrid*H.dyout ) ),
H.dyout )[:-1]
lons = np.arange( H.outlon0, ( H.outlon0 + ( H.numxgrid*H.dxout ) ),
H.dxout )[:-1]
data = data[:-1,:-1]
else:
lats = np.arange( H.outlat0, ( H.outlat0 + ( H.numygrid*H.dyout ) ),
H.dyout )
lons = np.arange( H.outlon0, ( H.outlon0 + ( H.numxgrid*H.dxout ) ),
H.dxout )
print data.shape
## transform to nx x ny regularly spaced native projection grid
if transform:
if m.projection not in ['cyl','merc','mill']:
dx = 2.*np.pi*m.rmajor/len(lons)
dy = 2.*np.pi*m.rminor/len(lats)
else:
dx = len(lons)
dy = len(lats)
nx = int((m.xmax-m.xmin)/dx)+1;
ny = int((m.ymax-m.ymin)/dy)+1
print nx
if nx is 1:
topodat = data
else:
topodat = m.transform_scalar(data,lons,lats,nx,ny)
else:
topodat = data
## Get the current axes, and properties for use later
pos = ax.get_position()
l, b, w, h = pos.bounds
## Set up the IMAGE
colmap = plt.get_cmap('gist_ncar')
im = m.imshow(topodat,cmap=colmap)
m.drawcoastlines()
return fig
class SuperDict(dict):
"""just so I can use . notation"""
def __getattr__(self, attr):
return self[attr]
def __setattr__(self, attr, value):
self[attr] = value
if __name__ == "__main__":
H = SuperDict()
H.name = 'working example'
H.outlat0 = -90
H.numygrid = 180
H.dyout = 1.
H.outlon0 = -179
H.numxgrid = 360
H.dxout = 1.0
H.data = np.random.rand(H.numygrid,H.numxgrid)
print H.data.shape
fig = plot_imshow_custom(H,transform=True)
plt.show()
print 'it worked'
try:
H.name = 'Not working example'
H.outlat0 = 40
H.numygrid = 100
H.dyout = 0.5
H.outlon0 = -179
H.numxgrid = 110
H.dxout = 0.5
H.data = np.random.rand(H.numygrid,H.numxgrid)
fig = plot_imshow_custom(H)
print 'huh?'
plt.show()
except:
print "As I said, it's not working..."
Jeff Whitaker wrote:
···
John [H2O] wrote:
I'm trying to 'automate' a few components within basemap. I have a pretty
complicated, and assuredly poorly written, set of functions that allow me
to
'dynamically' plot a grid of data (lon,lat).
Here is one section where I try to deal with transforming the data based
on
the projection. 'data' is a grid, often of size 720x360 or 720x180,
representing full globe or hemisphere at 0.5 degree resolution.
'outlon0',
outlat0', and 'd*out' are the llcrnr coordinates and step. 'transform' is
an
option, that is set to True by default:
1680 ## set up transformations for the data array
1681 if m.projection not in ['cyl','merc','mill']:
1682 lats = np.arange( outlat0, ( outlat0 + ( numygrid*dyout ) ),
dyout )[:-1]
1683 lons = np.arange( outlon0, ( outlon0 + ( numxgrid*dxout ) ),
dxout )[:-1]
1684 data = data[:-1,:-1]
1685 else:
1686 lats = np.arange( outlat0, ( outlat0 + ( numygrid*dyout ) ),
dyout )
1687 lons = np.arange( outlon0, ( outlon0 + ( numxgrid*dxout ) ),
dxout )
1688
1689 ## transform to nx x ny regularly spaced native projection grid
1690 if transform:
1691 dx = 2.*np.pi*m.rmajor/len(lons)
1692 nx = int((m.xmax-m.xmin)/dx)+1; ny =
int((m.ymax-m.ymin)/dx)+1
1693 if nx is 1:
1694 topodat = data
1695 else:
1696 topodat = m.transform_scalar(data,lons,lats,nx,ny)
1697 else:
1698 topodat = data
The problem is, when I use the approach with a 'cyl' grid, then
subsequently
try to draw the lsmask, I get a failure. Is this approach incorrect? I
had
to use the if nx is 1 statement because it was crashing with zero
division
error in some cases.
Thanks.
John: Please supply us with a self-contained example triggering the
error that we can run.
-Jeff
------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register
now!
http://p.sf.net/sfu/devconf
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
matplotlib-users List Signup and Options
--
View this message in context: http://www.nabble.com/basemap%2C-transform_scalar-question-tp25649437p25795985.html
Sent from the matplotlib - users mailing list archive at Nabble.com.