Jeff, Thank you again for the quick response.
Unfortunately, I need to output the data onto a regular lat/lon 0.5 degree grid for including in another project (this is the format we need, so I'll write the data out to a file). For plotting, however, I would be interested in seeing the 'easy' way to plot it. I've made a number of basemap plots, but honestly, I've made so many things complicated. An example with a projected grid dataset (when the available parameters for the projection are easily available -- as in this case) would be really helpful! Regards, john On Mon, Oct 25, 2010 at 10:24 PM, Jeff Whitaker <jsw...@fastmail.fm> wrote: > 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 > > -- 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
