Ortho projection fails at south pole in pcolormesh

[awarnock3]

trafficstars

I am plotting gridded satellite data using basemap 1.0.7 (I can't get 1.0.8 via pip yet). The grid is 0.25 degrees resolution, global coverage.

For polar views, I use pcolormesh. It works fine for the north polar view:

    bmap = Basemap(projection='ortho', lat_0=90, lon_0=270, resolution='c')
    show_map = bmap.pcolormesh(lons, lats, darray, vmin=vminval, vmax=vmaxval,
                               norm=norm, cmap=this_cmap, alpha=1.0, latlon=True)

but the corresponding map centered on the south pole plots garbage:

    bmap = Basemap(projection='ortho', lat_0=-90, lon_0=270, resolution='c')
    show_map = bmap.pcolormesh(lons, lats, darray, vmin=vminval, vmax=vmaxval,
                               norm=norm, cmap=this_cmap, alpha=1.0, latlon=True)

However, pcolor plots correctly at the south pole. This one is fine:

    bmap = Basemap(projection='ortho', lat_0=-90, lon_0=270, resolution='c')
    show_map = bmap.pcolor(lons, lats, darray, vmin=vminval, vmax=vmaxval,
                               norm=norm, cmap=this_cmap, alpha=1.0, latlon=True)

In this case, my gridded data array is pretty small and the time difference between pcolormesh and pcolor is negligible, but I'm guessing this is something that ought to be fixed, unless it already has been addressed in 1.0.8.

[pelson]

Thanks @awarnock3. Any chance you could provide a SSCCE example which allows us to see the problem without needing to have access to your data (i.e. use np.arange to produce coordinates)? That should speed up development time, as well as allowing us to suggest potential problems with your approach / potential options beyond pcolor.

Appreciate the time you've taken to raise this despite you having a workaround. Thanks again!

[awarnock3]

Probably. Let me see what I can come up with.

As a matter of fact, I do speak for my iPad.

On Oct 2, 2015, at 3:52 AM, Phil Elson [email protected] wrote:

Thanks @awarnock3. Any chance you could provide a SSCCE example which allows us to see the problem without needing to have access to your data (i.e. use np.arange to produce coordinates)? That should speed up development time, as well as allowing us to suggest potential problems with your approach / potential options beyond pcolor.

Appreciate the time you've taken to raise this despite you having a workaround. Thanks again!

— Reply to this email directly or view it on GitHub.

[awarnock3]

Here's one I threw together, along with example plots showing how pcolor() does the right thing at -90 but pcolormesh does not.

Thanks!

On 10/2/15 3:52 AM, Phil Elson wrote:

Thanks @awarnock3 https://github.com/awarnock3. Any chance you could provide a SSCCE example which allows us to see the problem without needing to have access to your data (i.e. use np.arange to produce coordinates)? That should speed up development time, as well as allowing us to suggest potential problems with your approach / potential options beyond pcolor.

Appreciate the time you've taken to raise this despite you having a workaround. Thanks again!

— Reply to this email directly or view it on GitHub https://github.com/matplotlib/basemap/issues/207#issuecomment-144948647.

Archie

-- Archie Warnock [email protected] -- A/WWW Enterprises www.awcubed.com -- As a matter of fact, I do speak for my employer.

[awarnock3]

I forgot to include the dependencies that generated the example.

basemap==1.0.7 matplotlib==1.4.3 numpy==1.9.3

On 10/2/15 3:52 AM, Phil Elson wrote:

Thanks @awarnock3 https://github.com/awarnock3. Any chance you could provide a SSCCE example which allows us to see the problem without needing to have access to your data (i.e. use np.arange to produce coordinates)? That should speed up development time, as well as allowing us to suggest potential problems with your approach / potential options beyond pcolor.

Appreciate the time you've taken to raise this despite you having a workaround. Thanks again!

— Reply to this email directly or view it on GitHub https://github.com/matplotlib/basemap/issues/207#issuecomment-144948647.

Archie

-- Archie Warnock [email protected] -- A/WWW Enterprises www.awcubed.com -- As a matter of fact, I do speak for my employer.

[guziy]

Hi @awarnock3:

I am curious what happens if you convert lats and lons to the projection coordinates yourself and plot without latlon=True ?

Thanks

[awarnock3]

Same thing. That was actually my first iteration. My naive guess is that whatever offset is done internally to avoid the singularity at the south pole isn't quite right. It works if you center at +90 but not at -90.

Additionally, simply changing "pcolormesh" to "pcolor" also fixes the issue. In my case, the grid size isn't too large and the performance isn't much different but I'll have other cases where the resolution is higher and the images will be much larger.

On 10/3/15 10:27 AM, Huziy Oleksandr (Sasha) wrote:

Hi @awarnock3 https://github.com/awarnock3:

I am curious what happens if you convert lats and lons to the projection coordinates yourself and plot without latlon=True ?

Thanks

— Reply to this email directly or view it on GitHub https://github.com/matplotlib/basemap/issues/207#issuecomment-145253259.

Archie

-- Archie Warnock [email protected] -- A/WWW Enterprises www.awcubed.com -- As a matter of fact, I do speak for my employer.

[guziy]

@awarnock3 Have you sent your example somewhere? I would like to play with it.

Cheers

[awarnock3]

Here you go. I've attached the example code and sample outputs.

On 10/4/15 1:07 PM, Huziy Oleksandr (Sasha) wrote:

@awarnock3 https://github.com/awarnock3 Have you sent your example somewhere? I would like to play with it.

Cheers

— Reply to this email directly or view it on GitHub https://github.com/matplotlib/basemap/issues/207#issuecomment-145366902.

Archie

-- Archie Warnock [email protected] -- A/WWW Enterprises www.awcubed.com -- As a matter of fact, I do speak for my employer.

[guziy]

Sorry @awarnock3, but I am still not able to lacate the script(s)...

[awarnock3]

Sorry - it appears that my attachments didn't get through. Here's the sample code and two images.

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt

def meshtest(darray):
    """This demonstrates the bug in pcolormesh.
    Define bmap with lat_0 = +90 and both plots are the same.  Define it
    with lat_0 = -90 and pcolor creates the expected result but pcolormesh
    does not.
    """
    plt.gcf()
    dpi = 100
    x_pixels = 360
    y_pixels = 180

    fig_x_inches = x_pixels / float(dpi)
    fig_y_inches = y_pixels / float(dpi)

    # Create the default figure at the right dimensions
    fig = plt.figure(figsize=(fig_x_inches, fig_y_inches), dpi=dpi)

    bmap = Basemap(projection='ortho', lat_0=-90, lon_0=-90, resolution='c')
    this_cmap = plt.get_cmap('PiYG')
    lon_array = np.arange(-180., 180., 1.)
    lat_array = np.arange(-90, 90., 1.)

    lons, lats = np.meshgrid(lon_array, lat_array)
    bmap.pcolormesh(lons, lats, darray, vmin=1.0, vmax=14.0, cmap=this_cmap,
                    alpha=1.0, latlon=True)
    fig.savefig('pcolormesh.png', dpi=dpi, transparent=1.0)

    bmap.pcolor(lons, lats, darray, vmin=1.0, vmax=14.0, cmap=this_cmap,
                alpha=1.0, latlon=True)
    fig.savefig('pcolor.png', dpi=dpi, transparent=1.0)

    plt.clf()  # clear the plot in case anything is cached
    plt.close('all')

    return

def make_data():
    """Make a highly colorful 2D image.  I'm sure there's a more pythonic
    way to do this but it works."""
    nrows = 180
    ncols = 360
    square_size = 20
    square_len = square_size**2
    nsquares = (nrows * ncols) / square_len
    darray = np.zeros((nrows, ncols))

    for square in range(nsquares):
        k = np.random.randint(16)
        nrow = square_size * int((square_size * square) / nrows)
        ncol = (square_size * square) % nrows
        darray[ncol:ncol + square_size, nrow:nrow + square_size] = k
    return darray

def main():
    darray = make_data()
    meshtest(darray)
    return

if __name__ == '__main__':
    main()

[WeatherGod]

Which version of matplotlib?

On Wed, Oct 7, 2015 at 6:33 AM, Archie Warnock [email protected] wrote:

Sorry - it appears that my attachments didn't get through. Here's the sample code and two images.

[image: pcolor] https://cloud.githubusercontent.com/assets/6179594/10335072/a49ad210-6cbc-11e5-84a6-9205a1b19bb0.png [image: pcolormesh] https://cloud.githubusercontent.com/assets/6179594/10335073/a4a0187e-6cbc-11e5-9a61-fe4e6cd71b1d.png

from mpl_toolkits.basemap import Basemap import numpy as np import matplotlib.pyplot as plt

def meshtest(darray): """This demonstrates the bug in pcolormesh. Define bmap with lat_0 = +90 and both plots are the same. Define it with lat_0 = -90 and pcolor creates the expected result but pcolormesh does not. """ plt.gcf() dpi = 100 x_pixels = 360 y_pixels = 180

fig_x_inches = x_pixels / float(dpi)
fig_y_inches = y_pixels / float(dpi)

# Create the default figure at the right dimensions
fig = plt.figure(figsize=(fig_x_inches, fig_y_inches), dpi=dpi)

bmap = Basemap(projection='ortho', lat_0=-90, lon_0=-90, resolution='c')
this_cmap = plt.get_cmap('PiYG')
lon_array = np.arange(-180., 180., 1.)
lat_array = np.arange(-90, 90., 1.)

lons, lats = np.meshgrid(lon_array, lat_array)
bmap.pcolormesh(lons, lats, darray, vmin=1.0, vmax=14.0, cmap=this_cmap,
                alpha=1.0, latlon=True)
fig.savefig('pcolormesh.png', dpi=dpi, transparent=1.0)

bmap.pcolor(lons, lats, darray, vmin=1.0, vmax=14.0, cmap=this_cmap,
            alpha=1.0, latlon=True)
fig.savefig('pcolor.png', dpi=dpi, transparent=1.0)

plt.clf()  # clear the plot in case anything is cached
plt.close('all')

return

def make_data(): """Make a highly colorful 2D image. I'm sure there's a more pythonic way to do this but it works.""" nrows = 180 ncols = 360 square_size = 20 square_len = square_size**2 nsquares = (nrows * ncols) / square_len darray = np.zeros((nrows, ncols))

for square in range(nsquares):
    k = np.random.randint(16)
    nrow = square_size * int((square_size * square) / nrows)
    ncol = (square_size * square) % nrows
    darray[ncol:ncol + square_size, nrow:nrow + square_size] = k
return darray

def main(): darray = make_data() meshtest(darray) return

if name == 'main': main()

— Reply to this email directly or view it on GitHub https://github.com/matplotlib/basemap/issues/207#issuecomment-146147578.

[awarnock3]

Matplotlib 1.4.3 Basemap 1.0.7 Numpy 1.9.3

As a matter of fact, I do speak for my iPhone.

On Oct 7, 2015, at 9:27 AM, Benjamin Root [email protected] wrote:

Which version of matplotlib?

On Wed, Oct 7, 2015 at 6:33 AM, Archie Warnock [email protected] wrote:

Sorry - it appears that my attachments didn't get through. Here's the sample code and two images.

[image: pcolor] https://cloud.githubusercontent.com/assets/6179594/10335072/a49ad210-6cbc-11e5-84a6-9205a1b19bb0.png [image: pcolormesh] https://cloud.githubusercontent.com/assets/6179594/10335073/a4a0187e-6cbc-11e5-9a61-fe4e6cd71b1d.png

from mpl_toolkits.basemap import Basemap import numpy as np import matplotlib.pyplot as plt

def meshtest(darray): """This demonstrates the bug in pcolormesh. Define bmap with lat_0 = +90 and both plots are the same. Define it with lat_0 = -90 and pcolor creates the expected result but pcolormesh does not. """ plt.gcf() dpi = 100 x_pixels = 360 y_pixels = 180

fig_x_inches = x_pixels / float(dpi) fig_y_inches = y_pixels / float(dpi)

Create the default figure at the right dimensions

fig = plt.figure(figsize=(fig_x_inches, fig_y_inches), dpi=dpi)

bmap = Basemap(projection='ortho', lat_0=-90, lon_0=-90, resolution='c') this_cmap = plt.get_cmap('PiYG') lon_array = np.arange(-180., 180., 1.) lat_array = np.arange(-90, 90., 1.)

lons, lats = np.meshgrid(lon_array, lat_array) bmap.pcolormesh(lons, lats, darray, vmin=1.0, vmax=14.0, cmap=this_cmap, alpha=1.0, latlon=True) fig.savefig('pcolormesh.png', dpi=dpi, transparent=1.0)

bmap.pcolor(lons, lats, darray, vmin=1.0, vmax=14.0, cmap=this_cmap, alpha=1.0, latlon=True) fig.savefig('pcolor.png', dpi=dpi, transparent=1.0)

plt.clf() # clear the plot in case anything is cached plt.close('all')

return

def make_data(): """Make a highly colorful 2D image. I'm sure there's a more pythonic way to do this but it works.""" nrows = 180 ncols = 360 square_size = 20 square_len = square_size**2 nsquares = (nrows * ncols) / square_len darray = np.zeros((nrows, ncols))

for square in range(nsquares): k = np.random.randint(16) nrow = square_size * int((square_size * square) / nrows) ncol = (square_size * square) % nrows darray[ncol:ncol + square_size, nrow:nrow + square_size] = k return darray

def main(): darray = make_data() meshtest(darray) return

if name == 'main': main()

— Reply to this email directly or view it on GitHub https://github.com/matplotlib/basemap/issues/207#issuecomment-146147578.

— Reply to this email directly or view it on GitHub.