Screenshot instructions:
Windows
Mac
Red Hat Linux
Ubuntu
Click URL instructions:
Rightclick on ad, choose "Copy Link", then paste here →
(This may not be possible with some types of ads)
From: Will Hewson <wh18@le...>  20100402 10:27:27
Attachments:
Message as HTML

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/Basemaporthographicprojectionplotdoesn%27trespectglobeboundarytp28117655p28117655.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
From: Will Hewson <wh18@le...>  20100402 10:27:28
Attachments:
Message as HTML

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/Basemaporthographicprojectionplotdoesn%27trespectglobeboundarytp28117654p28117654.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
From: Jeff Whitaker <jswhit@fa...>  20100402 11:39:07
Attachments:
Message as HTML

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/Basemaporthographicprojectionplotdoesn%27trespectglobeboundarytp28117654p28117654.html>; > Sent from the matplotlib  users mailing list archive > <http://old.nabble.com/matplotlibusersf2906.html>; at Nabble.com. > > >  > Download Intel® Parallel Studio Eval > Try the new software tools for yourself. Speed compiling, find bugs > proactively, and finetune applications for parallel performance. > See why Intel Parallel Studio got high marks during beta. > http://p.sf.net/sfu/intelswdev > > > _______________________________________________ > Matplotlibusers mailing list > Matplotlibusers@... > https://lists.sourceforge.net/lists/listinfo/matplotlibusers > 
From: Jeff Whitaker <jswhit@fa...>  20100402 12:16:56
Attachments:
text/plain

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.xmaxm.xmin)/40000.)+1; ny = int((m.ymaxm.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() 
From: Will Hewson <wh18@le...>  20100402 12:32:33

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.xmaxm.xmin)/40000.)+1; ny = int((m.ymaxm.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/Basemaporthographicprojectionplotdoesn%27trespectglobeboundarytp28117654p28118555.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
From: Jeff Whitaker <jswhit@fa...>  20100402 12:35:51

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.xmaxm.xmin)/40000.)+1; ny = int((m.ymaxm.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() >> >> >> > 
From: Jeff Whitaker <jswhit@fa...>  20100404 23:24:37

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 
From: Jeff Whitaker <jswhit@fa...>  20100405 12:11:31

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 finetune applications for parallel performance. >> See why Intel Parallel Studio got high marks during beta. >> http://p.sf.net/sfu/intelswdev >> _______________________________________________ >> Matplotlibusers mailing list >> Matplotlibusers@... >> https://lists.sourceforge.net/lists/listinfo/matplotlibusers >> >> >> > 
From: Will Hewson <wh18@le...>  20100405 13:25:28

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® Parallel Studio Eval >>> Try the new software tools for yourself. Speed compiling, find bugs >>> proactively, and finetune applications for parallel performance. >>> See why Intel Parallel Studio got high marks during beta. >>> http://p.sf.net/sfu/intelswdev >>> _______________________________________________ >>> Matplotlibusers mailing list >>> Matplotlibusers@... >>> https://lists.sourceforge.net/lists/listinfo/matplotlibusers >>> >>> >>> >> > > >  > Download Intel® Parallel Studio Eval > Try the new software tools for yourself. Speed compiling, find bugs > proactively, and finetune applications for parallel performance. > See why Intel Parallel Studio got high marks during beta. > http://p.sf.net/sfu/intelswdev > _______________________________________________ > Matplotlibusers mailing list > Matplotlibusers@... > https://lists.sourceforge.net/lists/listinfo/matplotlibusers > >  View this message in context: http://old.nabble.com/Basemaporthographicprojectionplotdoesn%27trespectglobeboundarytp28117654p28139978.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
From: Jeff Whitaker <jswhit@fa...>  20100405 14:30:21
Attachments:
hewson.py

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® Parallel Studio Eval >>>> Try the new software tools for yourself. Speed compiling, find bugs >>>> proactively, and finetune applications for parallel performance. >>>> See why Intel Parallel Studio got high marks during beta. >>>> http://p.sf.net/sfu/intelswdev >>>> _______________________________________________ >>>> Matplotlibusers mailing list >>>> Matplotlibusers@... >>>> https://lists.sourceforge.net/lists/listinfo/matplotlibusers >>>> >>>> >>>> >>>> >>> >>> >> >>  >> Download Intel® Parallel Studio Eval >> Try the new software tools for yourself. Speed compiling, find bugs >> proactively, and finetune applications for parallel performance. >> See why Intel Parallel Studio got high marks during beta. >> http://p.sf.net/sfu/intelswdev >> _______________________________________________ >> Matplotlibusers mailing list >> Matplotlibusers@... >> https://lists.sourceforge.net/lists/listinfo/matplotlibusers >> >> >> >  Jeffrey S. Whitaker Phone : (303)4976313 Meteorologist FAX : (303)4976449 NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@... 325 Broadway Office : Skaggs Research Cntr 1D113 Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg 
From: Will Hewson <wh18@le...>  20100405 15:19:41

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® Parallel Studio Eval >>>>> Try the new software tools for yourself. Speed compiling, find bugs >>>>> proactively, and finetune applications for parallel performance. >>>>> See why Intel Parallel Studio got high marks during beta. >>>>> http://p.sf.net/sfu/intelswdev >>>>> _______________________________________________ >>>>> Matplotlibusers mailing list >>>>> Matplotlibusers@... >>>>> https://lists.sourceforge.net/lists/listinfo/matplotlibusers >>>>> >>>>> >>>>> >>>>> >>>> >>>> >>> >>>  >>> Download Intel® Parallel Studio Eval >>> Try the new software tools for yourself. Speed compiling, find bugs >>> proactively, and finetune applications for parallel performance. >>> See why Intel Parallel Studio got high marks during beta. >>> http://p.sf.net/sfu/intelswdev >>> _______________________________________________ >>> Matplotlibusers mailing list >>> Matplotlibusers@... >>> https://lists.sourceforge.net/lists/listinfo/matplotlibusers >>> >>> >>> >> > > >  > Jeffrey S. Whitaker Phone : (303)4976313 > Meteorologist FAX : (303)4976449 > NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@... > 325 Broadway Office : Skaggs Research Cntr 1D113 > Boulder, CO, USA 803033328 Web : http://tinyurl.com/5telg > > > >  > Download Intel® Parallel Studio Eval > Try the new software tools for yourself. Speed compiling, find bugs > proactively, and finetune applications for parallel performance. > See why Intel Parallel Studio got high marks during beta. > http://p.sf.net/sfu/intelswdev > _______________________________________________ > Matplotlibusers mailing list > Matplotlibusers@... > https://lists.sourceforge.net/lists/listinfo/matplotlibusers > >  View this message in context: http://old.nabble.com/Basemaporthographicprojectionplotdoesn%27trespectglobeboundarytp28117654p28141052.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
From: Will Hewson <wh18@le...>  20100404 17:06:50

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.xmaxm.xmin)/40000.)+1; ny = int((m.ymaxm.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 finetune applications for parallel performance. > See why Intel Parallel Studio got high marks during beta. > http://p.sf.net/sfu/intelswdev > _______________________________________________ > Matplotlibusers mailing list > Matplotlibusers@... > https://lists.sourceforge.net/lists/listinfo/matplotlibusers > >  View this message in context: http://old.nabble.com/Basemaporthographicprojectionplotdoesn%27trespectglobeboundarytp28117654p28133659.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
From: Will Hewson <wh18@le...>  20100405 10:16:09

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 finetune applications for parallel performance. > See why Intel Parallel Studio got high marks during beta. > http://p.sf.net/sfu/intelswdev > _______________________________________________ > Matplotlibusers mailing list > Matplotlibusers@... > https://lists.sourceforge.net/lists/listinfo/matplotlibusers > >  View this message in context: http://old.nabble.com/Basemaporthographicprojectionplotdoesn%27trespectglobeboundarytp28117654p28138677.html Sent from the matplotlib  users mailing list archive at Nabble.com. 
Sign up for the SourceForge newsletter:
No, thanks