|
From: Jeff W. <js...@fa...> - 2010-04-02 11:39:07
|
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® 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
>
|