Closer, but still not quite right... not sure what I'm doing wrong?? 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 > -- 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
