On 10/25/10 2:28 PM, John wrote: > 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
John: If the data is on the native projection grid (xin, yin), just do x,y = np.meshgrid(xin,yin) # xin, yin are 1-d, x,y are 2-d (same shape as data) map.contourf(x,y,data,clevs) # clevs are contour levs. Making the plot this way and comparing after interpolating to a lat/lon grid should help you debug. -Jeff > 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 >> >> > > -- 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
