From: John [H2O] <washakie@gm...>  20090928 17:26:07

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.xmaxm.xmin)/dx)+1; ny = int((m.ymaxm.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.  View this message in context: http://www.nabble.com/basemap%2Ctransform_scalarquestiontp25649437p25649437.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
From: Jeff Whitaker <jswhit@fa...>  20090929 12:02:52

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.xmaxm.xmin)/dx)+1; ny = int((m.ymaxm.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 selfcontained example triggering the error that we can run. Jeff 
From: John [H2O] <washakie@gm...>  20091007 23:10:21

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.xmaxm.xmin)/dx)+1; ny = int((m.ymaxm.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.xmaxm.xmin)/dx)+1; ny = >> int((m.ymaxm.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 selfcontained 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 912, 2009. Register > now! > http://p.sf.net/sfu/devconf > _______________________________________________ > Matplotlibusers mailing list > Matplotlibusers@... > https://lists.sourceforge.net/lists/listinfo/matplotlibusers > >  View this message in context: http://www.nabble.com/basemap%2Ctransform_scalarquestiontp25649437p25795985.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
From: Jeff Whitaker <jswhit@fa...>  20091009 05:34:21

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.xmaxm.xmin)/dx)+1; > ny = int((m.ymaxm.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.xmaxm.xmin)/dx)+1; ny = >>> int((m.ymaxm.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 selfcontained 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 912, 2009. Register >> now! >> http://p.sf.net/sfu/devconf >> _______________________________________________ >> Matplotlibusers mailing list >> Matplotlibusers@... >> https://lists.sourceforge.net/lists/listinfo/matplotlibusers >> >> >> > > 
From: Eric Firing <efiring@ha...>  20091011 23:36:43

John [H2O] wrote: > > Jeff Whitaker wrote: >> >> 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 >> >> > > Jeff, > > I started using pcolormesh, but ran into the problem that I really required > a log scale color map for the plotting. As I recall, this could only be done > with imshow  but perhaps that has changed? Regardless, I've done some All the "mappables" like images, pcolortype plots (including pcolormesh), scatter, and collections can be created with a colormap and a norm. Using norm=LogNorm() gives you a log scale for color. See http://matplotlib.sourceforge.net/examples/pylab_examples/pcolor_log.html?highlight=pcolor_log > further testing. It seems I have the following problems with it: > > 1) It is much more 'coarse' in the coloring, imshow offers a finer/smoother > looking gradient. Imshow by default interpolates the colors; the pcolor* functions/methods do not. Eric > > 2) The bigger problem occurs when plotting in npstere, which I use quite > often. Maybe I have to add a wrap around, but as it is, I get a 'seam' at > the 180/180 meridian. > > 3) At the pole, there is a empty spot, creating a hole in my plots. This I > can't work with... :s > > Thanks again for your input. > > john > > 
From: John [H2O] <washakie@gm...>  20091011 21:11:03

Jeff Whitaker wrote: > > > 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 > > Jeff, I started using pcolormesh, but ran into the problem that I really required a log scale color map for the plotting. As I recall, this could only be done with imshow  but perhaps that has changed? Regardless, I've done some further testing. It seems I have the following problems with it: 1) It is much more 'coarse' in the coloring, imshow offers a finer/smoother looking gradient. 2) The bigger problem occurs when plotting in npstere, which I use quite often. Maybe I have to add a wrap around, but as it is, I get a 'seam' at the 180/180 meridian. 3) At the pole, there is a empty spot, creating a hole in my plots. This I can't work with... :s Thanks again for your input. john  View this message in context: http://www.nabble.com/basemap%2Ctransform_scalarquestiontp25649437p25847646.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
From: Jeff Whitaker <jswhit@fa...>  20091012 00:03:02

John [H2O] wrote: > Jeff Whitaker wrote: > >> 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 >> >> >> > > Jeff, > > I started using pcolormesh, but ran into the problem that I really required > a log scale color map for the plotting. As I recall, this could only be done > with imshow  but perhaps that has changed? John: The color mapping is the same with imshow or pcolormesh. See the recent thread on "logarithmic colormaps". > Regardless, I've done some > further testing. It seems I have the following problems with it: > > 1) It is much more 'coarse' in the coloring, imshow offers a finer/smoother > looking gradient. > Only way around that with pcolormesh is to interpolate to a finer grid. Imshow will interpolate the colors, pcolormesh will not. > 2) The bigger problem occurs when plotting in npstere, which I use quite > often. Maybe I have to add a wrap around, but as it is, I get a 'seam' at > the 180/180 meridian. > Just add a wraparound  that's easy with the addcyclic function. > 3) At the pole, there is a empty spot, creating a hole in my plots. This I > can't work with... :s > That's because you don't have a grid point at the pole. You could add one (by extrapolating the from the highest latitude). Jeff > Thanks again for your input. > > john > > > 