On 10/25/10 2:27 AM, John wrote: > 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 > > John: Basemap provides an 'interp' function that can be used to reproject data from one projection grid to another using bilinear interpolation.
http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html#mpl_toolkits.basemap.interp griddata is really for unstructured, scattered data. Since you have a regular grid, interp should be much faster. In your case, this should do it (untested). import numpy as np from mpl_toolkits.basemap import Basemap, interp # from http://nsidc.org/data/polar_stereo/ps_grids.html a = 6378.273e3 ec = 0.081816153 b = a*np.sqrt(1.-ec**2) map =\ Basemap(projection='stere',lat_0=90,lon_0=315,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 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.01,0.5) lats = np.arange(-90,90.01,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(datain, xin, yin, xout, yout,masked=True) -Jeff ------------------------------------------------------------------------------ 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
