John [H2O] wrote:
> 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.
>
John: I don't have time to look at your code right now, but let me just
make some general comments about plotting images on maps. If you want
to use imshow, the data your are plotting must coincide exactly with
your map plot area. So, for example if you want to plot a global
lat/lon grid on a north polar stereographic projection, you have to
interpolate to a rectangular grid in projection coordinates that fits in
the map region. However, in practice I find it's almost never worth
doing this. You can plot the data in the native coordinates on almost
any map projection region using pcolormesh or contourf, Just calculate
the x,y values of the of the data grid, and pass those values along with
the data to either one of those methods. Is there any particular reason
you want to use imshow, instead of pcolormesh or contourf?
-Jeff
> 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
>> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>>
>>
>>
>
>
------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) 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/devconference
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
