Jeff and others, I have a working solution and a non-working solution. I'm wondering whether someone can help me understand the problems in the latter. Please forgive in advance the long post, but all of my code is here.
First, you can get a dataset from : ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/ and the projection files here: http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats I've made a class for handling this data: <code> class NSIDC(object): """ Maybe a set of functions for NSIDC data """ def __init__(self,infile): self.infile = infile self.data = self.readseaice() def readseaice(self): """ reads the binary sea ice data and returns the header and the data see: http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html """ #use the BinaryFile class to access data from pflexpart import BinaryFile raw = BinaryFile(self.infile,'r') #start reading header values """ File Header Bytes Description 1-6 Missing data integer value 7-12 Number of columns in polar stereographic grid 13-18 Number of rows in polar stereographic grid 19-24 Unused/internal 25-30 Latitude enclosed by pointslar stereographic grid 31-36 Greenwich orientation of polar stereographicic grid 37-42 Unused/internal 43-48 J-coordinate of the grid intersection at the pole 49-54 I-coordinate of the grid intersection at the pole 55-60 Five-character instrument descriptor (SMMR, SSM/I) 61-66 Two descriptionriptors of two characters each that describe the data; for example, 07 cn = Nimbus-7 SMMR ice concentration 67-72 Starting Julian day of grid dayta 73-78 Starting hour of grid data (if available) 79-84 Starting minute of grid data (if available) 85-90 Ending Julian day of grid data 91-916 Ending hour of grid data (if available) 97-102 Ending minute of grid data (if available) 103-108 Year of grid data 109-114 Julian day of gridarea(xld data 115-120 Three-digit channel descriptor (000 for ice concentrationns) 121-126 Integer scaling factor 127-150 24-character file name 151-24 Unused3080-character image title 231-300 70-character data information (creation date, data source, etc.) """ hdr = raw.read(np.dtype('a1'),(300)) header = {} header['baddata'] = int(''.join(hdr[:6])) header['COLS'] = int(''.join(hdr[6:12])) header['ROWS'] = int(''.join(hdr[12:18])) header['lat'] = float(''.join(hdr[24:30])) header['lon0'] = float(''.join(hdr[30:36])) header['jcoord'] = float(''.join(hdr[42:48])) header['icoord'] = float(''.join(hdr[48:54])) header['instrument'] = ''.join(hdr[54:60]) header['descr'] = ''.join(hdr[60:66]) header['startjulday'] = int(''.join(hdr[66:72])) header['starthour'] = int(''.join(hdr[72:78])) header['startminute'] = int(''.join(hdr[78:84])) header['endjulday'] = int(''.join(hdr[84:90])) header['endhour'] = int(''.join(hdr[90:96])) header['endminute'] = int(''.join(hdr[96:102])) header['year'] = int(''.join(hdr[102:108])) header['julday'] = int(''.join(hdr[108:114])) header['chan'] = int(''.join(hdr[114:120])) header['scale'] = int(''.join(hdr[120:126])) header['filename'] = ''.join(hdr[126:150]) header['imagetitle'] = ''.join(hdr[150:230]) header['datainfo'] = ''.join(hdr[230:300]) #pdb.set_trace() seaiceconc = raw.read(np.uint8,(header['COLS'],header['ROWS'])) seaiceconc = seaiceconc.T return {'header':header,'data':seaiceconc} def conv2percentage(self): self.seaicepercentage = self.data['data']/2.5 def classify(self): """ classify the data into land, coast, missing, hole """ data = self.data['data'] self.header = self.data['header'] for a in [('land',254),('coast',253),('hole',251),('missing',255)]: zeros = np.zeros(data.shape) zeros[np.where(data==a[1])] = 1 exec('self.%s = zeros' % a[0]) #filter data data[data>250] = 0 self.ice = data def geocoordinate(self): """ use NSIDC grid files to assign lats/lons to grid. see: http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats """ try: ROWS = self.header['ROWS'] COLS = self.header['COLS'] except: raise AttributeError('object needs to have header, \ did you run self.classify?') datadir = '/flex_wrk/jfb/RESEARCH_ARCTIC/SEAICE/nsidc0081_ssmi_nrt_seaice' lonfile = os.path.join(datadir,'psn25lons_v2.dat') lons = np.fromfile(lonfile,dtype=np.dtype('i4'))/100000. self.lons = lons.reshape(ROWS,COLS) latfile = os.path.join(datadir,'psn25lats_v2.dat') lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000. self.lats = lats.reshape(ROWS,COLS) areafile = os.path.join(datadir,'psn25area_v2.dat') area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000. self.area = area.reshape(ROWS,COLS) def regrid2globe(self,dres=0.5): """ use parameters from http://nsidc.org/data/polar_stereo/ps_grids.html to regrid the data onto a lat/lon grid with degree resolution of dres """ a = 6378.273e3 ec = 0.081816153 b = a*np.sqrt(1.-ec**2) map = Basemap(projection='stere',lat_0=90,lon_0=-45,lat_ts=70,\ llcrnrlat=33.92,llcrnrlon=279.96,\ urcrnrlon=102.34,urcrnrlat=31.37,\ rsphere=(a,b)) # Basemap coordinate system starts with 0,0 at lower left corner nx = self.lons.shape[1] ny = self.lats.shape[0] xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points on the grid yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points on the grid # 0.5 degree grid lons = np.arange(-180,180,0.5) lats = np.arange(-90,90,0.5) lons, lats = np.meshgrid(lons,lats) xout,yout = map(lons, lats) # datain is the data on the nx,ny stereographic grid. # masked=True returns masked values for points outside projection grid dataout = interp(self.ice, xin, yin,\ xout, yout, masked=True) x,y = np.meshgrid(xin,yin) # xin, yin are 1-d, x,y are 2-d (same shape as data) map.contourf(x,y,np.flipud(self.ice),clevs=np.arange(0,250,1)) self.iceglobe = dataout self.m = map </code> You can use this as follows: import NSIDC d = NSIDC('nt_20080827_f17_nrt_n.bin') d.classify() d.geocoordinate() d.regrid2globe() d.m.drawcoastlines() This should produce a plot. QUESTION ONE: It seems possibly I'm slightly off (look at the great lakes). Any suggestions as to why? QUESTION TWO: Please, suggested improvements, code review or simplification is welcome. Now, on to the NON-WORKING example. I have a function (it's not too pretty, but it basically takes a bunch of pre-defined regions and plots a 'grid' of data to them. I've tried to make the necessary edits so you can use the code as is. data = d.iceglobe # using the d from the above example lons = np.arange(-180,180,0.5) lats = np.arange(-90,90,0.5) fig,m = plot_grid((lons,lats,data),region=None) QUESTION THREE: Notice that I'm quite sure the grid is upside down (you can see again the great lakes and Hudson bay outline in the top left). However, I've tried a variety of np.fliplr / np.flipud / data.T combinations, but I can't seem to get it projected properly. Any ideas here?? QUESTION FOUR: As above in Q2. <code> def plot_grid(D,region=None,dres=0.5, transform=True,figname=None,fillcontinents=False, points=False): """ plot an array over a basemap. The required argument "D" is either: * a tuple of (x,y,z) * a tuple with (lon,lat,grid) * a tuple with (grid,) Usage:: > FIG = plot_grid(D,**kwargs) =============== ======================================== keyword description =============== ======================================== dres resolution of the grid fillcontinents fills continents if True plotargs dictionary of plot arguments figname A figurename, if passed will save the figure with figname points set to True if passing a x,y,z matrix of points region A region from :func:`get_base1` =============== ======================================== """ print "length of D: %s" % len(D) if isinstance(D,np.ndarray): assert len(D.shape) == 2, "D grid must be 2d" #if len(D)==1: #assume full earth grid with dres print 'received grid of shape:', D.shape if D.shape[0] == 720: lons = np.arange(-180,180,dres) elif D.shape[0] == 721: lons = np.arange(-180,180.01,dres) if D.shape[1] == 360: lats = np.arange(-90,90,dres) elif D.shape[1] == 361: lats = np.arange(-90,90.01,dres) points = False z = D.T elif len(D)==3: points = True x = D[0] y = D[1] z = D[2] if len(z.shape) > 1: points = False lons = x lats = y #Set up a basemap ## CHANGED THIS HERE FOR THE CODEX ## ## Be sure to set region=None ## if isinstance(region,str): print "getting basemap with region: %s" % region fig,m = get_base1(region=region) else: fig = plt.figure() ax = fig.add_axes() a = 6378.273e3 ec = 0.081816153 b = a*np.sqrt(1.-ec**2) m = Basemap(projection='stere',lat_0=90,lon_0=-45,lat_ts=70,\ llcrnrlat=33.92,llcrnrlon=279.96,\ urcrnrlon=102.34,urcrnrlat=31.37,\ rsphere=(a,b)) m.drawcoastlines() if points == True: # Plot individual data points norm = colors.normalize(z.min(),z.max()) for i in range(len(y)): xpt,ypt = m(x[i],y[i]) cmap = cm.jet(norm(z[i])) #cmap = 'k' m.plot([xpt],[ypt],'.',color=cmap,markersize=2) #plt.plot(x[i],y[i],'.',color=cmap,markersize=2) if points == False: #transform Z data into projection #transform to nx x ny regularly spaced native projection grid dx = 2.*np.pi*m.rmajor/len(lons) nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1 if transform: # Need this if we choose lon,lat approach Zt,xx,yy = m.transform_scalar(z,lons,lats,nx,ny,returnxy=True) else: Zt = z print m.projection m.imshow(Zt) if fillcontinents: m.fillcontinents() plt.colorbar() #plt.imshow(Z) if figname != None: #plt.ylim([40,90]) #plt.title('data locations on mercator grid') plt.savefig(figname) else: plt.show() </code> ------------------------------------------------------------------------------ Nokia and AT&T present the 2010 Calling All Innovators-North America contest Create new apps & games for the Nokia N8 for consumers in U.S. and Canada $10 million total in prizes - $4M cash, 500 devices, nearly $6M in marketing Develop with Nokia Qt SDK, Web Runtime, or Java and Publish to Ovi Store http://p.sf.net/sfu/nokia-dev2dev _______________________________________________ Matplotlib-users mailing list Matplotlib-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/matplotlib-users
