matplotlib-users

[Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Will H. <wh...@le...>]

Hi forum/ mailing list,

When I plot in the orthographic projection I'm getting the large artefact
shown below extending away from the north east of the globe.

I'm not finding the same problem when plotting in a full globe projection so
I'm presuming the problem is with the way I'm projecting everything rather
than my data itself.

I've included my plotting code below, if anyone is able to spot some glaring
omissions/ errors I'd be most grateful (I've been using python/ matplotlib
for only a couple of weeks now!).


http://old.nabble.com/file/p28117655/binploterr.png 


#!/usr/local/bin/python2.6

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.basemap import Basemap

import sys, glob


#input must be 3 col file of lons lats and data

#bins input values into half degree grid, ignores negative values


plts = glob.glob('*.plt')

x = np.arange(-180, 180, 0.5); y = np.arange(-90, 90, 0.5)

grid_lon, grid_lat = np.meshgrid(x,y) #regularly spaced 2D grid

n_vals = np.zeros((360,720)) #mean divisor

dat = np.zeros((360,720)) #2D grid of zeros 


for pt in plts:

    

    in_file = pt

    data = np.loadtxt(in_file, comments = ';')

    fname = in_file.split('.')[0]


    lon = data[:,0] #original 1D list

    lat = data[:,1] #original 1D list

    slcol = data[:,2] #z data


    lon = (np.around(lon*2))/2 #round to nearest .0 or 0.5

    lat = (np.around(lat*2))/2 #round to nearest .0 or 0.5


    ##keep the below between files


    j=0


    for i in slcol:

        if lon[j] < 0:

            grid_lon_ind = 360+(lon[j]*2)

            grid_lat_ind = 180+(lat[j]*2)

        else:

            grid_lon_ind = 360-(lon[j]*2)

            grid_lat_ind = 180+(lat[j]*2)

            

        if i > 0:

            dat[grid_lat_ind, grid_lon_ind] += i #add i'th value

            n_vals[grid_lat_ind, grid_lon_ind] += 1 #increase cell counter
by 1 for each extra value

        j+=1


dat = np.nan_to_num(dat/n_vals)


#create map object

fig = plt.figure()

m = Basemap(projection='ortho', lon_0=lon[(len(lon)/2)], lat_0=0,
resolution='l', area_thresh=10000.)

#m = Basemap(projection='moll',lon_0=0,resolution='c', area_thresh=10000.)


X,Y = m(grid_lon, grid_lat)


#pass all 2d arrays to pcolor

im = m.pcolormesh(X,Y,dat)


#add coastlines, globe boundary and colourbar

m.drawcoastlines()

m.drawmapboundary()

m.drawparallels(np.arange(-90, 90,30))

m.drawmeridians(np.arange(-180,180,30))


fig.colorbar(im)

plt.title('CH20 and ting')

plt.savefig('binplot.png')



Thanks for your help,


Will.
-- 
View this message in context: http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117655p28117655.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

[Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Will H. <wh...@le...>]

Hi forum/ mailing list,

When I plot in the orthographic projection I'm getting the large artefact
shown below extending away from the north east of the globe.

I'm not finding the same problem when plotting in a full globe projection so
I'm presuming the problem is with the way I'm projecting everything rather
than my data itself.

I've included my plotting code below, if anyone is able to spot some glaring
omissions/ errors I'd be most grateful (I've been using python/ matplotlib
for only a couple of weeks now!).


http://old.nabble.com/file/p28117654/binploterr.png 


#!/usr/local/bin/python2.6

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.basemap import Basemap

import sys, glob


#input must be 3 col file of lons lats and data

#bins input values into half degree grid, ignores negative values


plts = glob.glob('*.plt')

x = np.arange(-180, 180, 0.5); y = np.arange(-90, 90, 0.5)

grid_lon, grid_lat = np.meshgrid(x,y) #regularly spaced 2D grid

n_vals = np.zeros((360,720)) #mean divisor

dat = np.zeros((360,720)) #2D grid of zeros 


for pt in plts:

    

    in_file = pt

    data = np.loadtxt(in_file, comments = ';')

    fname = in_file.split('.')[0]


    lon = data[:,0] #original 1D list

    lat = data[:,1] #original 1D list

    slcol = data[:,2] #z data


    lon = (np.around(lon*2))/2 #round to nearest .0 or 0.5

    lat = (np.around(lat*2))/2 #round to nearest .0 or 0.5


    ##keep the below between files


    j=0


    for i in slcol:

        if lon[j] < 0:

            grid_lon_ind = 360+(lon[j]*2)

            grid_lat_ind = 180+(lat[j]*2)

        else:

            grid_lon_ind = 360-(lon[j]*2)

            grid_lat_ind = 180+(lat[j]*2)

            

        if i > 0:

            dat[grid_lat_ind, grid_lon_ind] += i #add i'th value

            n_vals[grid_lat_ind, grid_lon_ind] += 1 #increase cell counter
by 1 for each extra value

        j+=1


dat = np.nan_to_num(dat/n_vals)


#create map object

fig = plt.figure()

m = Basemap(projection='ortho', lon_0=lon[(len(lon)/2)], lat_0=0,
resolution='l', area_thresh=10000.)

#m = Basemap(projection='moll',lon_0=0,resolution='c', area_thresh=10000.)


X,Y = m(grid_lon, grid_lat)


#pass all 2d arrays to pcolor

im = m.pcolormesh(X,Y,dat)


#add coastlines, globe boundary and colourbar

m.drawcoastlines()

m.drawmapboundary()

m.drawparallels(np.arange(-90, 90,30))

m.drawmeridians(np.arange(-180,180,30))


fig.colorbar(im)

plt.title('CH20 and ting')

plt.savefig('binplot.png')



Thanks for your help,


Will.
-- 
View this message in context: http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117654p28117654.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Jeff W. <js...@fa...>]

On 4/2/10 4:27 AM, Will Hewson wrote:
> Hi forum/ mailing list, When I plot in the orthographic projection I'm 
> getting the large artefact shown below extending away from the north 
> east of the globe. I'm not finding the same problem when plotting in a 
> full globe projection so I'm presuming the problem is with the way I'm 
> projecting everything rather than my data itself. I've included my 
> plotting code below, if anyone is able to spot some glaring omissions/ 
> errors I'd be most grateful (I've been using python/ matplotlib for 
> only a couple of weeks now!).
>

Will:   You'll have to provide the data so we can actually run the script.

-Jeff
> #!/usr/local/bin/python2.6
> import numpy as np
> import matplotlib.pyplot as plt
> from mpl_toolkits.basemap import Basemap
> import sys, glob
>
> #input must be 3 col file of lons lats and data
> #bins input values into half degree grid, ignores negative values
>
> plts = glob.glob('*.plt')
> x = np.arange(-180, 180, 0.5); y = np.arange(-90, 90, 0.5)
> grid_lon, grid_lat = np.meshgrid(x,y) #regularly spaced 2D grid
> n_vals = np.zeros((360,720)) #mean divisor
> dat = np.zeros((360,720)) #2D grid of zeros
>
> for pt in plts:
>
> in_file = pt
> data = np.loadtxt(in_file, comments = ';')
> fname = in_file.split('.')[0]
>
> lon = data[:,0] #original 1D list
> lat = data[:,1] #original 1D list
> slcol = data[:,2] #z data
>
> lon = (np.around(lon*2))/2 #round to nearest .0 or 0.5
> lat = (np.around(lat*2))/2 #round to nearest .0 or 0.5
>
> ##keep the below between files
>
> j=0
>
> for i in slcol:
> if lon[j] < 0:
> grid_lon_ind = 360+(lon[j]*2)
> grid_lat_ind = 180+(lat[j]*2)
> else:
> grid_lon_ind = 360-(lon[j]*2)
> grid_lat_ind = 180+(lat[j]*2)
>
> if i > 0:
> dat[grid_lat_ind, grid_lon_ind] += i #add i'th value
> n_vals[grid_lat_ind, grid_lon_ind] += 1 #increase cell counter by 1 
> for each extra value
> j+=1
>
> dat = np.nan_to_num(dat/n_vals)
>
> #create map object
> fig = plt.figure()
> m = Basemap(projection='ortho', lon_0=lon[(len(lon)/2)], lat_0=0, 
> resolution='l', area_thresh=10000.)
> #m = Basemap(projection='moll',lon_0=0,resolution='c', area_thresh=10000.)
>
> X,Y = m(grid_lon, grid_lat)
>
> #pass all 2d arrays to pcolor
> im = m.pcolormesh(X,Y,dat)
>
> #add coastlines, globe boundary and colourbar
> m.drawcoastlines()
> m.drawmapboundary()
> m.drawparallels(np.arange(-90, 90,30))
> m.drawmeridians(np.arange(-180,180,30))
>
> fig.colorbar(im)
> plt.title('CH20 and ting')
> plt.savefig('binplot.png')
>
> Thanks for your help,
> Will.
> ------------------------------------------------------------------------
> View this message in context: Basemap/ orthographic projection plot 
> doesn't respect globe boundary 
> <http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117654p28117654.html>
> Sent from the matplotlib - users mailing list archive 
> <http://old.nabble.com/matplotlib---users-f2906.html> at Nabble.com.
>
>
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
>
>
> _______________________________________________
> Matplotlib-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>    

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Jeff W. <js...@fa...>]

from mpl_toolkits.basemap import Basemap, shiftgrid
import numpy as np
import matplotlib.pyplot as plt
# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
topoin = np.loadtxt('etopo20data.gz')
lons = np.loadtxt('etopo20lons.gz')
lats = np.loadtxt('etopo20lats.gz')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lons = shiftgrid(180.,topoin,lons,start=False)
m = Basemap(projection='ortho',lon_0=-105,lat_0=40,resolution='l')
# transform to nx x ny regularly spaced native projection grid
nx = int((m.xmax-m.xmin)/40000.)+1; ny = int((m.ymax-m.ymin)/40000.)+1
topodat,x,y =\
m.transform_scalar(topoin,lons,lats,nx,ny,returnxy=True,masked=True,order=1)
# create the figure.
fig=plt.figure(figsize=(8,8))
im = m.pcolormesh(x,y,topodat,cmap=plt.cm.jet)
m.drawcoastlines()
m.drawparallels(np.arange(0.,80,20.))
m.drawmeridians(np.arange(10.,360.,30.))
m.drawmapboundary()
plt.show()

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Will H. <wh...@le...>]

This is great Jeff, thanks for the help - I'll give it a try over the weekend
(it's bank holiday here in the UK!) and get back to you, if I'm still having
trouble I'll stick up the plotting data too... thanks again.

Will



Jeff Whitaker wrote:
> 
> On 4/2/10 4:27 AM, Will Hewson wrote:
>> Hi forum/ mailing list, When I plot in the orthographic projection I'm 
>> getting the large artefact shown below extending away from the north 
>> east of the globe. I'm not finding the same problem when plotting in a 
>> full globe projection so I'm presuming the problem is with the way I'm 
>> projecting everything rather than my data itself. I've included my 
>> plotting code below, if anyone is able to spot some glaring omissions/ 
>> errors I'd be most grateful (I've been using python/ matplotlib for 
>> only a couple of weeks now!).
> Will:  I think what's happening is that pcolormesh is having trouble 
> dealing with the higher curvlinear grid, which becomes nearly 
> pathological near the horizon of the projection.  If you take a look at 
> the test.py file in the basemap examples directory, you'll see an 
> example orthographic plot that solves this problem by first 
> interpolating the data to a regular grid in projection coordinates (with 
> values over the plot horizon masked).  The example uses imshow, but 
> pcolormesh works as well.  A standalone version of the example using 
> pcolormesah  is attached, which uses data files in the basemap examples 
> directory.
> 
> -Jeff
> 
> from mpl_toolkits.basemap import Basemap, shiftgrid
> import numpy as np
> import matplotlib.pyplot as plt
> # read in topo data (on a regular lat/lon grid)
> # longitudes go from 20 to 380.
> topoin = np.loadtxt('etopo20data.gz')
> lons = np.loadtxt('etopo20lons.gz')
> lats = np.loadtxt('etopo20lats.gz')
> # shift data so lons go from -180 to 180 instead of 20 to 380.
> topoin,lons = shiftgrid(180.,topoin,lons,start=False)
> m = Basemap(projection='ortho',lon_0=-105,lat_0=40,resolution='l')
> # transform to nx x ny regularly spaced native projection grid
> nx = int((m.xmax-m.xmin)/40000.)+1; ny = int((m.ymax-m.ymin)/40000.)+1
> topodat,x,y =\
> m.transform_scalar(topoin,lons,lats,nx,ny,returnxy=True,masked=True,order=1)
> # create the figure.
> fig=plt.figure(figsize=(8,8))
> im = m.pcolormesh(x,y,topodat,cmap=plt.cm.jet)
> m.drawcoastlines()
> m.drawparallels(np.arange(0.,80,20.))
> m.drawmeridians(np.arange(10.,360.,30.))
> m.drawmapboundary()
> plt.show()
> 
> 

-- 
View this message in context: http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117654p28118555.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Jeff W. <js...@fa...>]

On 4/2/10 6:32 AM, Will Hewson wrote:
> This is great Jeff, thanks for the help - I'll give it a try over the weekend
> (it's bank holiday here in the UK!) and get back to you, if I'm still having
> trouble I'll stick up the plotting data too... thanks again.
>
> Will
>    

Will:  I forgot to mention that contourf will work on your data without 
having to interpolate to projection coordinates.

-Jeff
>
>
> Jeff Whitaker wrote:
>    
>> On 4/2/10 4:27 AM, Will Hewson wrote:
>>      
>>> Hi forum/ mailing list, When I plot in the orthographic projection I'm
>>> getting the large artefact shown below extending away from the north
>>> east of the globe. I'm not finding the same problem when plotting in a
>>> full globe projection so I'm presuming the problem is with the way I'm
>>> projecting everything rather than my data itself. I've included my
>>> plotting code below, if anyone is able to spot some glaring omissions/
>>> errors I'd be most grateful (I've been using python/ matplotlib for
>>> only a couple of weeks now!).
>>>        
>> Will:  I think what's happening is that pcolormesh is having trouble
>> dealing with the higher curvlinear grid, which becomes nearly
>> pathological near the horizon of the projection.  If you take a look at
>> the test.py file in the basemap examples directory, you'll see an
>> example orthographic plot that solves this problem by first
>> interpolating the data to a regular grid in projection coordinates (with
>> values over the plot horizon masked).  The example uses imshow, but
>> pcolormesh works as well.  A standalone version of the example using
>> pcolormesah  is attached, which uses data files in the basemap examples
>> directory.
>>
>> -Jeff
>>
>> from mpl_toolkits.basemap import Basemap, shiftgrid
>> import numpy as np
>> import matplotlib.pyplot as plt
>> # read in topo data (on a regular lat/lon grid)
>> # longitudes go from 20 to 380.
>> topoin = np.loadtxt('etopo20data.gz')
>> lons = np.loadtxt('etopo20lons.gz')
>> lats = np.loadtxt('etopo20lats.gz')
>> # shift data so lons go from -180 to 180 instead of 20 to 380.
>> topoin,lons = shiftgrid(180.,topoin,lons,start=False)
>> m = Basemap(projection='ortho',lon_0=-105,lat_0=40,resolution='l')
>> # transform to nx x ny regularly spaced native projection grid
>> nx = int((m.xmax-m.xmin)/40000.)+1; ny = int((m.ymax-m.ymin)/40000.)+1
>> topodat,x,y =\
>> m.transform_scalar(topoin,lons,lats,nx,ny,returnxy=True,masked=True,order=1)
>> # create the figure.
>> fig=plt.figure(figsize=(8,8))
>> im = m.pcolormesh(x,y,topodat,cmap=plt.cm.jet)
>> m.drawcoastlines()
>> m.drawparallels(np.arange(0.,80,20.))
>> m.drawmeridians(np.arange(10.,360.,30.))
>> m.drawmapboundary()
>> plt.show()
>>
>>
>>      
>    

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Jeff W. <js...@fa...>]

On 4/4/10 11:06 AM, Will Hewson wrote:
> Hi again Jeff et al...
>
> I've had a play around with the extra few lines of code - on paper this
> seems like it should solve the problems I'm experiencing. However, an
> error's being thrown up by the transform scalar function, as my lons and
> lats won't necessarily be increasing. The data I'm plotting is satellite
> data and so at the beginning and end of the orbit file lats go over the pole
> from 90 to -90, with a similar problem for the lons - whereby the data is
> taken across the satellite track. I've thought about sorting the data before
> passing it to transform_scalar but I'm always going to be left with the
> problem in either lats or lons.
>
> I've uploaded the file I'm currently working with this time. It's three
> columns of lons, lats and z values.
>
> Once again, many thanks for your help.
>
> Will.
>
> http://old.nabble.com/file/p28133659/test.plt test.plt
>    

Will:  Is it a regular lat/lon grid or a satellite swath?  If it's the 
latter, you can't use my solution.

-Jeff

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Jeff W. <js...@fa...>]

On 4/5/10 4:16 AM, Will Hewson wrote:
> Hey Jeff,
>
> It's somewhere between the two - the original satellite swath is converted
> to a regular 0.5 degree grid by truncating, binning, and averaging each
> point's lons and lats over the top of a 720 x 360 np.zeros array. the
> plotting still works fine for non ortho/ hemispherical projections, and I've
> no big problem with using global projections for the time being. Thanks for
> your help in the meantime anyway.
>
> Cheers,
>
>
> Will.
>    

Will:  If it's a regular 0.5 degree lat/lon grid, it should work in 
transform_scalar.  However, I don't see how to read the data in your 
test.plt file into a regular 360x720 grid.  It seems to only contain the 
points in the swath with nonzero values.

-Jeff
>
>
> Jeff Whitaker wrote:
>    
>> On 4/4/10 11:06 AM, Will Hewson wrote:
>>      
>>> Hi again Jeff et al...
>>>
>>> I've had a play around with the extra few lines of code - on paper this
>>> seems like it should solve the problems I'm experiencing. However, an
>>> error's being thrown up by the transform scalar function, as my lons and
>>> lats won't necessarily be increasing. The data I'm plotting is satellite
>>> data and so at the beginning and end of the orbit file lats go over the
>>> pole
>>> from 90 to -90, with a similar problem for the lons - whereby the data is
>>> taken across the satellite track. I've thought about sorting the data
>>> before
>>> passing it to transform_scalar but I'm always going to be left with the
>>> problem in either lats or lons.
>>>
>>> I've uploaded the file I'm currently working with this time. It's three
>>> columns of lons, lats and z values.
>>>
>>> Once again, many thanks for your help.
>>>
>>> Will.
>>>
>>> http://old.nabble.com/file/p28133659/test.plt test.plt
>>>
>>>        
>> Will:  Is it a regular lat/lon grid or a satellite swath?  If it's the
>> latter, you can't use my solution.
>>
>> -Jeff
>>
>>
>> ------------------------------------------------------------------------------
>> Download Intel® Parallel Studio Eval
>> Try the new software tools for yourself. Speed compiling, find bugs
>> proactively, and fine-tune applications for parallel performance.
>> See why Intel Parallel Studio got high marks during beta.
>> http://p.sf.net/sfu/intel-sw-dev
>> _______________________________________________
>> Matplotlib-users mailing list
>> Mat...@li...
>> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>>
>>
>>      
>    

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Will H. <wh...@le...>]

I should perhaps of explained my code (included in top post) a little better,
the values in my attached file aren't on a regular grid to start with, I do
a little bit of juggling as follows to get them into a regular grid:

I'm firstly setting up my 2D grid of 0.5 degree lat lons, followed by
identically sized grids of zeros for the data bin, and mean divisors:

x = np.arange(-180, 180, 0.5); y = np.arange(-90, 90, 0.5)
grid_lon, grid_lat = np.meshgrid(x,y) #regularly spaced 2D grid
n_vals = np.zeros((360,720)) #mean divisor
dat = np.zeros((360,720)) #2D grid of zeros 

I'm then taking my input data (e.g. the .plt file attached), and rounding
the lat lons to the nearest 0 or 0.5:

lon = (np.around(lon*2))/2 #round to nearest .0 or 0.5
lat = (np.around(lat*2))/2 #round to nearest .0 or 0.5

Then for each row in my input file where Z is greater than 0, I'm adding the
n'th Z value to its corresponding position in the dat zeros array, and
keeping a count of how many values are going into each cell in the mean
divisor array:

j=0
for i in slcol:
  if lon[j] < 0:
    grid_lon_ind = 360+(lon[j]*2)
    grid_lat_ind = 180+(lat[j]*2)
  else:
    grid_lon_ind = 360-(lon[j]*2)
    grid_lat_ind = 180+(lat[j]*2)
  if i > 0:
    dat[grid_lat_ind, grid_lon_ind] += i #add i'th value
    n_vals[grid_lat_ind, grid_lon_ind] += 1 #increase cell counter by 1 for
each extra value
  j+=1

Finally the new dat array is divided by the mean divisor array to give me my
mean Z values:

dat = np.nan_to_num(dat/n_vals)

I've done it this way as opposed to interpolating *properly* in order to
(for instance) stop the values bleeding away from the edges of the satellite
swath.

Cheers,

Will.


Jeff Whitaker wrote:
> 
> On 4/5/10 4:16 AM, Will Hewson wrote:
>> Hey Jeff,
>>
>> It's somewhere between the two - the original satellite swath is
>> converted
>> to a regular 0.5 degree grid by truncating, binning, and averaging each
>> point's lons and lats over the top of a 720 x 360 np.zeros array. the
>> plotting still works fine for non ortho/ hemispherical projections, and
>> I've
>> no big problem with using global projections for the time being. Thanks
>> for
>> your help in the meantime anyway.
>>
>> Cheers,
>>
>>
>> Will.
>>    
> 
> Will:  If it's a regular 0.5 degree lat/lon grid, it should work in 
> transform_scalar.  However, I don't see how to read the data in your 
> test.plt file into a regular 360x720 grid.  It seems to only contain the 
> points in the swath with nonzero values.
> 
> -Jeff
>>
>>
>> Jeff Whitaker wrote:
>>    
>>> On 4/4/10 11:06 AM, Will Hewson wrote:
>>>      
>>>> Hi again Jeff et al...
>>>>
>>>> I've had a play around with the extra few lines of code - on paper this
>>>> seems like it should solve the problems I'm experiencing. However, an
>>>> error's being thrown up by the transform scalar function, as my lons
>>>> and
>>>> lats won't necessarily be increasing. The data I'm plotting is
>>>> satellite
>>>> data and so at the beginning and end of the orbit file lats go over the
>>>> pole
>>>> from 90 to -90, with a similar problem for the lons - whereby the data
>>>> is
>>>> taken across the satellite track. I've thought about sorting the data
>>>> before
>>>> passing it to transform_scalar but I'm always going to be left with the
>>>> problem in either lats or lons.
>>>>
>>>> I've uploaded the file I'm currently working with this time. It's three
>>>> columns of lons, lats and z values.
>>>>
>>>> Once again, many thanks for your help.
>>>>
>>>> Will.
>>>>
>>>> http://old.nabble.com/file/p28133659/test.plt test.plt
>>>>
>>>>        
>>> Will:  Is it a regular lat/lon grid or a satellite swath?  If it's the
>>> latter, you can't use my solution.
>>>
>>> -Jeff
>>>
>>>
>>> ------------------------------------------------------------------------------
>>> Download Intel&#174; Parallel Studio Eval
>>> Try the new software tools for yourself. Speed compiling, find bugs
>>> proactively, and fine-tune applications for parallel performance.
>>> See why Intel Parallel Studio got high marks during beta.
>>> http://p.sf.net/sfu/intel-sw-dev
>>> _______________________________________________
>>> Matplotlib-users mailing list
>>> Mat...@li...
>>> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>>>
>>>
>>>      
>>    
> 
> 
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Matplotlib-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
> 
> 

-- 
View this message in context: http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117654p28139978.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Jeff W. <js...@fa...>]

On 4/5/10 7:25 AM, Will Hewson wrote:
> I should perhaps of explained my code (included in top post) a little better,
> the values in my attached file aren't on a regular grid to start with, I do
> a little bit of juggling as follows to get them into a regular grid:
>
> I'm firstly setting up my 2D grid of 0.5 degree lat lons, followed by
> identically sized grids of zeros for the data bin, and mean divisors:
>
> x = np.arange(-180, 180, 0.5); y = np.arange(-90, 90, 0.5)
> grid_lon, grid_lat = np.meshgrid(x,y) #regularly spaced 2D grid
> n_vals = np.zeros((360,720)) #mean divisor
> dat = np.zeros((360,720)) #2D grid of zeros
>
> I'm then taking my input data (e.g. the .plt file attached), and rounding
> the lat lons to the nearest 0 or 0.5:
>
> lon = (np.around(lon*2))/2 #round to nearest .0 or 0.5
> lat = (np.around(lat*2))/2 #round to nearest .0 or 0.5
>
> Then for each row in my input file where Z is greater than 0, I'm adding the
> n'th Z value to its corresponding position in the dat zeros array, and
> keeping a count of how many values are going into each cell in the mean
> divisor array:
>
> j=0
> for i in slcol:
>    if lon[j]<  0:
>      grid_lon_ind = 360+(lon[j]*2)
>      grid_lat_ind = 180+(lat[j]*2)
>    else:
>      grid_lon_ind = 360-(lon[j]*2)
>      grid_lat_ind = 180+(lat[j]*2)
>    if i>  0:
>      dat[grid_lat_ind, grid_lon_ind] += i #add i'th value
>      n_vals[grid_lat_ind, grid_lon_ind] += 1 #increase cell counter by 1 for
> each extra value
>    j+=1
>
> Finally the new dat array is divided by the mean divisor array to give me my
> mean Z values:
>
> dat = np.nan_to_num(dat/n_vals)
>
> I've done it this way as opposed to interpolating *properly* in order to
> (for instance) stop the values bleeding away from the edges of the satellite
> swath.
>
> Cheers,
>
> Will.
>    

Will:  I made some slight modifications to your original script and it 
works fine with the ortho projection using either contourf on the 
original lat/lon grid or pcolormesh on the interpolated map projection 
grid.

-Jeff
>
> Jeff Whitaker wrote:
>    
>> On 4/5/10 4:16 AM, Will Hewson wrote:
>>      
>>> Hey Jeff,
>>>
>>> It's somewhere between the two - the original satellite swath is
>>> converted
>>> to a regular 0.5 degree grid by truncating, binning, and averaging each
>>> point's lons and lats over the top of a 720 x 360 np.zeros array. the
>>> plotting still works fine for non ortho/ hemispherical projections, and
>>> I've
>>> no big problem with using global projections for the time being. Thanks
>>> for
>>> your help in the meantime anyway.
>>>
>>> Cheers,
>>>
>>>
>>> Will.
>>>
>>>        
>> Will:  If it's a regular 0.5 degree lat/lon grid, it should work in
>> transform_scalar.  However, I don't see how to read the data in your
>> test.plt file into a regular 360x720 grid.  It seems to only contain the
>> points in the swath with nonzero values.
>>
>> -Jeff
>>      
>>>
>>> Jeff Whitaker wrote:
>>>
>>>        
>>>> On 4/4/10 11:06 AM, Will Hewson wrote:
>>>>
>>>>          
>>>>> Hi again Jeff et al...
>>>>>
>>>>> I've had a play around with the extra few lines of code - on paper this
>>>>> seems like it should solve the problems I'm experiencing. However, an
>>>>> error's being thrown up by the transform scalar function, as my lons
>>>>> and
>>>>> lats won't necessarily be increasing. The data I'm plotting is
>>>>> satellite
>>>>> data and so at the beginning and end of the orbit file lats go over the
>>>>> pole
>>>>> from 90 to -90, with a similar problem for the lons - whereby the data
>>>>> is
>>>>> taken across the satellite track. I've thought about sorting the data
>>>>> before
>>>>> passing it to transform_scalar but I'm always going to be left with the
>>>>> problem in either lats or lons.
>>>>>
>>>>> I've uploaded the file I'm currently working with this time. It's three
>>>>> columns of lons, lats and z values.
>>>>>
>>>>> Once again, many thanks for your help.
>>>>>
>>>>> Will.
>>>>>
>>>>> http://old.nabble.com/file/p28133659/test.plt test.plt
>>>>>
>>>>>
>>>>>            
>>>> Will:  Is it a regular lat/lon grid or a satellite swath?  If it's the
>>>> latter, you can't use my solution.
>>>>
>>>> -Jeff
>>>>
>>>>
>>>> ------------------------------------------------------------------------------
>>>> Download Intel&#174; Parallel Studio Eval
>>>> Try the new software tools for yourself. Speed compiling, find bugs
>>>> proactively, and fine-tune applications for parallel performance.
>>>> See why Intel Parallel Studio got high marks during beta.
>>>> http://p.sf.net/sfu/intel-sw-dev
>>>> _______________________________________________
>>>> Matplotlib-users mailing list
>>>> Mat...@li...
>>>> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>>>>
>>>>
>>>>
>>>>          
>>>
>>>        
>>
>> ------------------------------------------------------------------------------
>> Download Intel&#174; Parallel Studio Eval
>> Try the new software tools for yourself. Speed compiling, find bugs
>> proactively, and fine-tune applications for parallel performance.
>> See why Intel Parallel Studio got high marks during beta.
>> http://p.sf.net/sfu/intel-sw-dev
>> _______________________________________________
>> Matplotlib-users mailing list
>> Mat...@li...
>> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>>
>>
>>      
>    


-- 
Jeffrey S. Whitaker         Phone  : (303)497-6313
Meteorologist               FAX    : (303)497-6449
NOAA/OAR/PSD  R/PSD1        Email  : Jef...@no...
325 Broadway                Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web    : http://tinyurl.com/5telg

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Will H. <wh...@le...>]

Jeff, this is great, works fine - many thanks for all your help over the last
few days, it really is appreciated. I'm trying to build the case within my
office for switching over to Basemap from IDL, ironing out niggles like this
is really useful in this respect.

All the best,


Will.


Jeff Whitaker wrote:
> 
> On 4/5/10 7:25 AM, Will Hewson wrote:
>> I should perhaps of explained my code (included in top post) a little
>> better,
>> the values in my attached file aren't on a regular grid to start with, I
>> do
>> a little bit of juggling as follows to get them into a regular grid:
>>
>> I'm firstly setting up my 2D grid of 0.5 degree lat lons, followed by
>> identically sized grids of zeros for the data bin, and mean divisors:
>>
>> x = np.arange(-180, 180, 0.5); y = np.arange(-90, 90, 0.5)
>> grid_lon, grid_lat = np.meshgrid(x,y) #regularly spaced 2D grid
>> n_vals = np.zeros((360,720)) #mean divisor
>> dat = np.zeros((360,720)) #2D grid of zeros
>>
>> I'm then taking my input data (e.g. the .plt file attached), and rounding
>> the lat lons to the nearest 0 or 0.5:
>>
>> lon = (np.around(lon*2))/2 #round to nearest .0 or 0.5
>> lat = (np.around(lat*2))/2 #round to nearest .0 or 0.5
>>
>> Then for each row in my input file where Z is greater than 0, I'm adding
>> the
>> n'th Z value to its corresponding position in the dat zeros array, and
>> keeping a count of how many values are going into each cell in the mean
>> divisor array:
>>
>> j=0
>> for i in slcol:
>>    if lon[j]<  0:
>>      grid_lon_ind = 360+(lon[j]*2)
>>      grid_lat_ind = 180+(lat[j]*2)
>>    else:
>>      grid_lon_ind = 360-(lon[j]*2)
>>      grid_lat_ind = 180+(lat[j]*2)
>>    if i>  0:
>>      dat[grid_lat_ind, grid_lon_ind] += i #add i'th value
>>      n_vals[grid_lat_ind, grid_lon_ind] += 1 #increase cell counter by 1
>> for
>> each extra value
>>    j+=1
>>
>> Finally the new dat array is divided by the mean divisor array to give me
>> my
>> mean Z values:
>>
>> dat = np.nan_to_num(dat/n_vals)
>>
>> I've done it this way as opposed to interpolating *properly* in order to
>> (for instance) stop the values bleeding away from the edges of the
>> satellite
>> swath.
>>
>> Cheers,
>>
>> Will.
>>    
> 
> Will:  I made some slight modifications to your original script and it 
> works fine with the ortho projection using either contourf on the 
> original lat/lon grid or pcolormesh on the interpolated map projection 
> grid.
> 
> -Jeff
>>
>> Jeff Whitaker wrote:
>>    
>>> On 4/5/10 4:16 AM, Will Hewson wrote:
>>>      
>>>> Hey Jeff,
>>>>
>>>> It's somewhere between the two - the original satellite swath is
>>>> converted
>>>> to a regular 0.5 degree grid by truncating, binning, and averaging each
>>>> point's lons and lats over the top of a 720 x 360 np.zeros array. the
>>>> plotting still works fine for non ortho/ hemispherical projections, and
>>>> I've
>>>> no big problem with using global projections for the time being. Thanks
>>>> for
>>>> your help in the meantime anyway.
>>>>
>>>> Cheers,
>>>>
>>>>
>>>> Will.
>>>>
>>>>        
>>> Will:  If it's a regular 0.5 degree lat/lon grid, it should work in
>>> transform_scalar.  However, I don't see how to read the data in your
>>> test.plt file into a regular 360x720 grid.  It seems to only contain the
>>> points in the swath with nonzero values.
>>>
>>> -Jeff
>>>      
>>>>
>>>> Jeff Whitaker wrote:
>>>>
>>>>        
>>>>> On 4/4/10 11:06 AM, Will Hewson wrote:
>>>>>
>>>>>          
>>>>>> Hi again Jeff et al...
>>>>>>
>>>>>> I've had a play around with the extra few lines of code - on paper
>>>>>> this
>>>>>> seems like it should solve the problems I'm experiencing. However, an
>>>>>> error's being thrown up by the transform scalar function, as my lons
>>>>>> and
>>>>>> lats won't necessarily be increasing. The data I'm plotting is
>>>>>> satellite
>>>>>> data and so at the beginning and end of the orbit file lats go over
>>>>>> the
>>>>>> pole
>>>>>> from 90 to -90, with a similar problem for the lons - whereby the
>>>>>> data
>>>>>> is
>>>>>> taken across the satellite track. I've thought about sorting the data
>>>>>> before
>>>>>> passing it to transform_scalar but I'm always going to be left with
>>>>>> the
>>>>>> problem in either lats or lons.
>>>>>>
>>>>>> I've uploaded the file I'm currently working with this time. It's
>>>>>> three
>>>>>> columns of lons, lats and z values.
>>>>>>
>>>>>> Once again, many thanks for your help.
>>>>>>
>>>>>> Will.
>>>>>>
>>>>>> http://old.nabble.com/file/p28133659/test.plt test.plt
>>>>>>
>>>>>>
>>>>>>            
>>>>> Will:  Is it a regular lat/lon grid or a satellite swath?  If it's the
>>>>> latter, you can't use my solution.
>>>>>
>>>>> -Jeff
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>> Download Intel&#174; Parallel Studio Eval
>>>>> Try the new software tools for yourself. Speed compiling, find bugs
>>>>> proactively, and fine-tune applications for parallel performance.
>>>>> See why Intel Parallel Studio got high marks during beta.
>>>>> http://p.sf.net/sfu/intel-sw-dev
>>>>> _______________________________________________
>>>>> Matplotlib-users mailing list
>>>>> Mat...@li...
>>>>> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>>>>>
>>>>>
>>>>>
>>>>>          
>>>>
>>>>        
>>>
>>> ------------------------------------------------------------------------------
>>> Download Intel&#174; Parallel Studio Eval
>>> Try the new software tools for yourself. Speed compiling, find bugs
>>> proactively, and fine-tune applications for parallel performance.
>>> See why Intel Parallel Studio got high marks during beta.
>>> http://p.sf.net/sfu/intel-sw-dev
>>> _______________________________________________
>>> Matplotlib-users mailing list
>>> Mat...@li...
>>> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>>>
>>>
>>>      
>>    
> 
> 
> -- 
> Jeffrey S. Whitaker         Phone  : (303)497-6313
> Meteorologist               FAX    : (303)497-6449
> NOAA/OAR/PSD  R/PSD1        Email  : Jef...@no...
> 325 Broadway                Office : Skaggs Research Cntr 1D-113
> Boulder, CO, USA 80303-3328 Web    : http://tinyurl.com/5telg
> 
> 
>  
> ------------------------------------------------------------------------------
> Download Intel&#174; Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Matplotlib-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
> 
> 

-- 
View this message in context: http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117654p28141052.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Will H. <wh...@le...>]

Hi again Jeff et al...

I've had a play around with the extra few lines of code - on paper this
seems like it should solve the problems I'm experiencing. However, an
error's being thrown up by the transform scalar function, as my lons and
lats won't necessarily be increasing. The data I'm plotting is satellite
data and so at the beginning and end of the orbit file lats go over the pole
from 90 to -90, with a similar problem for the lons - whereby the data is
taken across the satellite track. I've thought about sorting the data before
passing it to transform_scalar but I'm always going to be left with the
problem in either lats or lons.

I've uploaded the file I'm currently working with this time. It's three
columns of lons, lats and z values.

Once again, many thanks for your help.

Will.

http://old.nabble.com/file/p28133659/test.plt test.plt 


Jeff Whitaker wrote:
> 
> On 4/2/10 6:32 AM, Will Hewson wrote:
>> This is great Jeff, thanks for the help - I'll give it a try over the
>> weekend
>> (it's bank holiday here in the UK!) and get back to you, if I'm still
>> having
>> trouble I'll stick up the plotting data too... thanks again.
>>
>> Will
>>    
> 
> Will:  I forgot to mention that contourf will work on your data without 
> having to interpolate to projection coordinates.
> 
> -Jeff
>>
>>
>> Jeff Whitaker wrote:
>>    
>>> On 4/2/10 4:27 AM, Will Hewson wrote:
>>>      
>>>> Hi forum/ mailing list, When I plot in the orthographic projection I'm
>>>> getting the large artefact shown below extending away from the north
>>>> east of the globe. I'm not finding the same problem when plotting in a
>>>> full globe projection so I'm presuming the problem is with the way I'm
>>>> projecting everything rather than my data itself. I've included my
>>>> plotting code below, if anyone is able to spot some glaring omissions/
>>>> errors I'd be most grateful (I've been using python/ matplotlib for
>>>> only a couple of weeks now!).
>>>>        
>>> Will:  I think what's happening is that pcolormesh is having trouble
>>> dealing with the higher curvlinear grid, which becomes nearly
>>> pathological near the horizon of the projection.  If you take a look at
>>> the test.py file in the basemap examples directory, you'll see an
>>> example orthographic plot that solves this problem by first
>>> interpolating the data to a regular grid in projection coordinates (with
>>> values over the plot horizon masked).  The example uses imshow, but
>>> pcolormesh works as well.  A standalone version of the example using
>>> pcolormesah  is attached, which uses data files in the basemap examples
>>> directory.
>>>
>>> -Jeff
>>>
>>> from mpl_toolkits.basemap import Basemap, shiftgrid
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> # read in topo data (on a regular lat/lon grid)
>>> # longitudes go from 20 to 380.
>>> topoin = np.loadtxt('etopo20data.gz')
>>> lons = np.loadtxt('etopo20lons.gz')
>>> lats = np.loadtxt('etopo20lats.gz')
>>> # shift data so lons go from -180 to 180 instead of 20 to 380.
>>> topoin,lons = shiftgrid(180.,topoin,lons,start=False)
>>> m = Basemap(projection='ortho',lon_0=-105,lat_0=40,resolution='l')
>>> # transform to nx x ny regularly spaced native projection grid
>>> nx = int((m.xmax-m.xmin)/40000.)+1; ny = int((m.ymax-m.ymin)/40000.)+1
>>> topodat,x,y =\
>>> m.transform_scalar(topoin,lons,lats,nx,ny,returnxy=True,masked=True,order=1)
>>> # create the figure.
>>> fig=plt.figure(figsize=(8,8))
>>> im = m.pcolormesh(x,y,topodat,cmap=plt.cm.jet)
>>> m.drawcoastlines()
>>> m.drawparallels(np.arange(0.,80,20.))
>>> m.drawmeridians(np.arange(10.,360.,30.))
>>> m.drawmapboundary()
>>> plt.show()
>>>
>>>
>>>      
>>    
> 
> 
> ------------------------------------------------------------------------------
> Download Intel® Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Matplotlib-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
> 
> 

-- 
View this message in context: http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117654p28133659.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

[Will H. <wh...@le...>]

Hey Jeff,

It's somewhere between the two - the original satellite swath is converted
to a regular 0.5 degree grid by truncating, binning, and averaging each
point's lons and lats over the top of a 720 x 360 np.zeros array. the
plotting still works fine for non ortho/ hemispherical projections, and I've
no big problem with using global projections for the time being. Thanks for
your help in the meantime anyway.

Cheers,


Will.



Jeff Whitaker wrote:
> 
> On 4/4/10 11:06 AM, Will Hewson wrote:
>> Hi again Jeff et al...
>>
>> I've had a play around with the extra few lines of code - on paper this
>> seems like it should solve the problems I'm experiencing. However, an
>> error's being thrown up by the transform scalar function, as my lons and
>> lats won't necessarily be increasing. The data I'm plotting is satellite
>> data and so at the beginning and end of the orbit file lats go over the
>> pole
>> from 90 to -90, with a similar problem for the lons - whereby the data is
>> taken across the satellite track. I've thought about sorting the data
>> before
>> passing it to transform_scalar but I'm always going to be left with the
>> problem in either lats or lons.
>>
>> I've uploaded the file I'm currently working with this time. It's three
>> columns of lons, lats and z values.
>>
>> Once again, many thanks for your help.
>>
>> Will.
>>
>> http://old.nabble.com/file/p28133659/test.plt test.plt
>>    
> 
> Will:  Is it a regular lat/lon grid or a satellite swath?  If it's the 
> latter, you can't use my solution.
> 
> -Jeff
> 
> 
> ------------------------------------------------------------------------------
> Download Intel® Parallel Studio Eval
> Try the new software tools for yourself. Speed compiling, find bugs
> proactively, and fine-tune applications for parallel performance.
> See why Intel Parallel Studio got high marks during beta.
> http://p.sf.net/sfu/intel-sw-dev
> _______________________________________________
> Matplotlib-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
> 
> 

-- 
View this message in context: http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117654p28138677.html
Sent from the matplotlib - users mailing list archive at Nabble.com.