Apologies, I see I didn't need to change the xin, yin variables in the interp function. I have it working now, but I still don't quite have the plotting correct... be back with a report. -john
On Mon, Oct 25, 2010 at 9:27 PM, John <washa...@gmail.com> wrote: > Jeff, thanks for the answer. Unfortunately I have a problem due to the > 'polar' nature of the grid, the xin, yin values are not increasing. I > tried both passing lat and lon grids or flattened vectors of the > original data, but neither works. You can see my method here, which is > a new method of my NSIDC class: > > 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=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 > nx = self.lons.shape[0] > 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.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(self.ice.flatten(), self.lons.flatten(), > self.lats.flatten(),\ > xout, yout,masked=True) > > self.regridded = dataout > > Thank you, > john > > On Mon, Oct 25, 2010 at 1:51 PM, Jeff Whitaker <jsw...@fastmail.fm> wrote: >> 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 >> >> >> >> > > > > -- > Configuration > `````````````````````````` > Plone 2.5.3-final, > CMF-1.6.4, > Zope (Zope 2.9.7-final, python 2.4.4, linux2), > Python 2.6 > PIL 1.1.6 > Mailman 2.1.9 > Postfix 2.4.5 > Procmail v3.22 2001/09/10 > -- Configuration `````````````````````````` Plone 2.5.3-final, CMF-1.6.4, Zope (Zope 2.9.7-final, python 2.4.4, linux2), Python 2.6 PIL 1.1.6 Mailman 2.1.9 Postfix 2.4.5 Procmail v3.22 2001/09/10 ------------------------------------------------------------------------------ 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
