Hello, I'm trying to take a npstereographic grid of points and reproject it so that I can save out a file in regular 0.5 degree lat lon grid coordinates. The description of the grid points in the current projection is here: http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
I've written the following class for handling the data: 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 francesc 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'])) 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 = '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(COLS,ROWS) latfile = os.path.join(datadir,'psn25lats_v2.dat') lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000. self.lats = lats.reshape(COLS,ROWS) areafile = os.path.join(datadir,'psn25area_v2.dat') area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000. self.area = area.reshape(COLS,ROWS) Once I have the data in python, I've done some plotting with some weird results, I'm clearly not doing something correctly. I'd like to know the best way to do this, and I suspect it would require GDAL. Here's what I'm trying to do: from NSIDC import NSIDC import numpy as np from matplotlib import mlab, pyplot as plt #load the data d = jfb.NSIDC('nsidc0081_ssmi_nrt_seaice/nt_20080827_f17_nrt_n.bin') d.classify() d.geocoordinate() x = d.lons.ravel(); y = d.lats.ravel(); z = d.ice.ravel() #create a regular lat/lon grid 0.5 degrees dres=0.5 reg_lon = np.arange(-180,180,dres,'f') reg_lat=np.arange(-90,90,dres,'f') #regrid the data into the regular grid Z = mlab.griddata(x,y,z,reg_lon,reg_lat) My result is confusing: plotting the data as imshow, demonstrates I'm loading it okay: plt.imshow(d.ice) yields: http://picasaweb.google.com/washakie/Temp#5531888474069952818 however, if I try to plot the reprojected grid on to a map, I get something strange: http://picasaweb.google.com/washakie/Temp#5531888480458500386 But if I use the same map object to plot my original data using the scatter function: m.scatter(x,y,c=z) I also get a strange result, so possibly I'm not reading the grid in properly, but it seems strange, since all the 'points' are where they are supposed to be, but I would not expect this color pattern: http://picasaweb.google.com/washakie/Temp#5531895787146118914 Is anyone willing to write a quick tutorial or script on how this would be achieved properly in gdal or basemap? Thanks, john ------------------------------------------------------------------------------ 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