On 10/25/10 2:13 PM, John wrote: > Closer, but still not quite right... not sure what I'm doing wrong??
John: Since I don't know what the plot should look like, it's hard to say. Perhaps the data is just transposed? Let me back up a bit and ask why you want to interpolate to a lat/lon grid? If it's just to make a plot, you can plot with Basemap on the original projection grid quite easily. -Jeff > http://picasaweb.google.com/lh/photo/6Dnylo0BcjX0A-wdmczmNg?feat=directlink > > Any ideas?? > > -john > > On Mon, Oct 25, 2010 at 9:53 PM, John<washa...@gmail.com> wrote: >> 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 >> > > -- Jeffrey S. Whitaker Phone : (303)497-6313 Meteorologist FAX : (303)497-6449 NOAA/OAR/PSD R/PSD1 Email : jeffrey.s.whita...@noaa.gov 325 Broadway Office : Skaggs Research Cntr 1D-113 Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg ------------------------------------------------------------------------------ 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
