Basemap/ orthographic projection plot doesn't respect globe boundary

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

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

Sent from the matplotlib - users mailing list archive at Nabble.com.

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

···

http://p.sf.net/sfu/intel-sw-dev


Matplotlib-users@lists.sourceforge.nethttps://lists.sourceforge.net/lists/listinfo/matplotlib-users

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

plotmap_ortho.py (930 Bytes)

···

http://p.sf.net/sfu/intel-sw-dev


Matplotlib-users@lists.sourceforge.nethttps://lists.sourceforge.net/lists/listinfo/matplotlib-users

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't-respect-globe-boundary-tp28117654p28118555.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

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

···

On 4/2/10 6:32 AM, Will Hewson wrote:

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()

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


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

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

-Jeff

···

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

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


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

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

···

On 4/5/10 4:16 AM, Will Hewson wrote:

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

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


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

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

hewson.py (1.94 KB)

···

On 4/5/10 7:25 AM, Will Hewson wrote:

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
Matplotlib-users@lists.sourceforge.net
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
Matplotlib-users@lists.sourceforge.net
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 : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

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
Matplotlib-users@lists.sourceforge.net
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
Matplotlib-users@lists.sourceforge.net
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 : Jeffrey.S.Whitaker@...259...
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
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users


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