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