Nicolas, I've attached a script I used to load the files, metadata and coordinates.
You owe me a donut. David On Fri, Feb 26, 2010 at 8:11 AM, Dag Sverre Seljebotn <da...@student.matnat.uio.no> wrote: > Nicolas wrote: >> Hello >> >> a VERY specific question, but I thought someone in the list might have >> used this dataset before and could give me some advice >> >> I am trying to read the daily Antarctic Sea Ice Concentrations from >> Nimbus-7 SMMR and DMSP SSM/I Passive Microwave Data, as provided by >> the national snow and ice data center (http://nsidc.org), more >> specifically, these are the files as downloaded from their ftp site >> (ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/final-gsfc/south/daily) >> >> they are provided in binary files (e.g. nt_19980101_f13_v01_s.bin for >> the 1st of Jan. 1998) >> >> the metadata information >> (http://nsidc.org/cgi-bin/get_metadata.pl?id=nsidc-0051) gives the >> following information (for the polar stereographic projection): >> """ >> Data are scaled and stored as one-byte integers in flat binary arrays >> >> geographical coordinates >> N: 90° S: 30.98° E: 180° W: -180° >> >> Latitude Resolution: 25 km >> Longitude Resolution: 25 km >> >> Distribution Size: 105212 bytes per southern file >> """ >> >> I am unfamiliar with non self-documented files (used to hdf and netcdf >> !) and struggle to make sense of how to read these files and plot the >> corresponding maps, I've tried using the array module and the >> fromstring function >> >> file=open(filename,'rb') >> a=array.array('B',file.read()) >> var=numpy.fromstring(a,dtype=np.int) >> > Try > > numpy.fromfile(f, dtype=np.int8) > > or uint8, depending on whether the data is signed or not. What you did > is also correct except that the final dtype must be np.int8. > > (You can always do var = var.astype(np.int) *afterwards* to convert to a > bigger integer type.) > > Dag Sverre > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion >
""" Load data from SMMR and SSM/I ICE CONCENTRATION datasets. See http://nsidc.org/data/polar_stereo/tools.html The first 300 bytes contain the header. The data is stored as 1 byte integers. The number of columns is stored in the header[6:11] and the number of rows in header[12:17] :author: David Huard <david.hu...@gmail.com> :date: September 2009 Notes ===== ltln_25n.msk : Latitude and longitude lines, ie crosshairs. Cells containing land, coast, the hole at the pole or missing data are indicated as follows: land : 254 coast : 253 hole : 251 missing1 : 255 Resolution ---------- 25km : (448, 304) """ import numpy as np import struct, os from datetime import datetime, timedelta PATH = '/home/data/seaice/' def load_meta(filename): """Return the meta data about the file.""" meta = {} f = open(filename, 'r') header = f.read(126) meta['filename'] = f.read(24).strip(' \x00') meta['title'] = f.read(80).strip(' \x00') meta['info'] = f.read(70).strip(' \x00') h = [header[i:i+5] for i in range(0,122,6)] meta['missing'] = int(h[0]) meta['ncols'] = int(h[1]) meta['nrows'] = int(h[2]) meta['latitude'] = float(h[4]) meta['orientation'] = float(h[5]) # With respect to Greenwich meta['pole_j'] = int(float(h[7])) meta['pole_i'] = int(float(h[8])) meta['instrument'] = h[9] meta['data descriptors'] = h[10] meta['julianday_start'] = h[11] meta['hour_start'] = h[12] meta['minute_start'] = h[13] meta['julianday_end'] = h[14] meta['hour_end'] = int(h[15]) meta['minute_end'] = int(h[16]) meta['year'] = int(h[17]) meta['julianday'] = int(h[18]) meta['channel'] = h[19] meta['scaling_factor'] = int(h[20]) meta['packing'] = float(h[20]) try: meta['date_start'] = datetime(int(h[17]), 1, 1, int(h[12]), int(h[13])) + timedelta(days=h[11]) except: pass return meta def read(filename): meta = load_meta(filename) s = (meta['nrows'], meta['ncols']) x = np.memmap(filename, dtype=np.uint8, mode='r', shape=s, offset=300) return x def coords25(): """Return the longitude and latitude of the stereographic grid for 25 km resolution. """ s = (448, 304) lon = np.memmap(os.path.join(PATH, 'tools', 'psn25lons_v2.dat'), dtype='<i4', mode='r', shape=s) lat = np.memmap(os.path.join(PATH, 'tools', 'psn25lats_v2.dat'), dtype='<i4', mode='r', shape=s) return lon/100000., lat/100000. def iceplot(a): """Plot an ice concentration array of the CAA. """ import matplotlib as mpl import matplotlib.pyplot as plt sl = slice(130, 350), slice(40, 250), a = a[sl] land = a == 254 coast = a == 253 hole = a == 251 missing = a == 255 # Plot the coast fig = plt.figure() ax = fig.add_subplot(111) ax.imshow(np.ma.array(a, mask=np.logical_not(coast)), mpl.cm.gray) # Plot the ice concentration im = ax.imshow(np.ma.array(a/250., mask=a>250).round(2)) ax.set_xticklabels([0]) ax.set_yticklabels([0]) cb = plt.colorbar(mappable=im, ticks=np.linspace(0,1,11)) cb.set_clim(vmax=1) cb.set_label('Sea ice concentration') return fig
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
