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
------------------------------------------------------------------------------
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
