basemap, transform_scalar question

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.

···

--
View this message in context: http://www.nabble.com/basemap%2C-transform_scalar-question-tp25649437p25649437.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

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

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
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

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

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

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%2C-transform_scalar-question-tp25649437p25847646.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

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, pcolor-type 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

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