reprojecting a grid / mlab or basemap or gdal??

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

···

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

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 <jswhit@...146...> 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\):
    &quot;&quot;&quot; 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
    &quot;&quot;&quot;
    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\[&#39;COLS&#39;\],header\[&#39;ROWS&#39;\]\)\)

    return \{&#39;header&#39;:header,&#39;data&#39;:seaiceconc\}

def conv2percentage\(self\):
    self\.seaicepercentage = self\.data\[&#39;data&#39;\]/2\.5

def classify\(self\):
    &quot;&quot;&quot; classify the data into land, coast, missing, hole &quot;&quot;&quot;
    data = self\.data\[&#39;data&#39;\]
    self\.header = self\.data\[&#39;header&#39;\]

    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&gt;250\] = 0
    self\.ice = data

def geocoordinate\(self\):
    &quot;&quot;&quot; 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\[&#39;ROWS&#39;\]
        COLS = self\.header\[&#39;COLS&#39;\]
    except:
        raise AttributeError\(&#39;object needs to have header, \\
                did you run self\.classify?&#39;\)

    datadir = &#39;nsidc0081\_ssmi\_nrt\_seaice&#39;

    lonfile = os\.path\.join\(datadir,&#39;psn25lons\_v2\.dat&#39;\)
    lons = np\.fromfile\(lonfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/100000\.
    self\.lons = lons\.reshape\(COLS,ROWS\)

    latfile = os\.path\.join\(datadir,&#39;psn25lats\_v2\.dat&#39;\)
    lats = np\.fromfile\(latfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/100000\.
    self\.lats = lats\.reshape\(COLS,ROWS\)

    areafile = os\.path\.join\(datadir,&#39;psn25area\_v2\.dat&#39;\)
    area = np\.fromfile\(latfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/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

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 <washakie@...287...> 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=&#39;stere&#39;,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 <jswhit@...146...> 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\):
    &quot;&quot;&quot; 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
    &quot;&quot;&quot;
    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\[&#39;COLS&#39;\],header\[&#39;ROWS&#39;\]\)\)

    return \{&#39;header&#39;:header,&#39;data&#39;:seaiceconc\}

def conv2percentage\(self\):
    self\.seaicepercentage = self\.data\[&#39;data&#39;\]/2\.5

def classify\(self\):
    &quot;&quot;&quot; classify the data into land, coast, missing, hole &quot;&quot;&quot;
    data = self\.data\[&#39;data&#39;\]
    self\.header = self\.data\[&#39;header&#39;\]

    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&gt;250\] = 0
    self\.ice = data

def geocoordinate\(self\):
    &quot;&quot;&quot; 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\[&#39;ROWS&#39;\]
        COLS = self\.header\[&#39;COLS&#39;\]
    except:
        raise AttributeError\(&#39;object needs to have header, \\
                did you run self\.classify?&#39;\)

    datadir = &#39;nsidc0081\_ssmi\_nrt\_seaice&#39;

    lonfile = os\.path\.join\(datadir,&#39;psn25lons\_v2\.dat&#39;\)
    lons = np\.fromfile\(lonfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/100000\.
    self\.lons = lons\.reshape\(COLS,ROWS\)

    latfile = os\.path\.join\(datadir,&#39;psn25lats\_v2\.dat&#39;\)
    lats = np\.fromfile\(latfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/100000\.
    self\.lats = lats\.reshape\(COLS,ROWS\)

    areafile = os\.path\.join\(datadir,&#39;psn25area\_v2\.dat&#39;\)
    area = np\.fromfile\(latfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/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

</details>

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 <washakie@...287...> 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 <washakie@...287...> 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=&#39;stere&#39;,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 <jswhit@...146...> 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\):
    &quot;&quot;&quot; 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
    &quot;&quot;&quot;
    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\[&#39;COLS&#39;\],header\[&#39;ROWS&#39;\]\)\)

    return \{&#39;header&#39;:header,&#39;data&#39;:seaiceconc\}

def conv2percentage\(self\):
    self\.seaicepercentage = self\.data\[&#39;data&#39;\]/2\.5

def classify\(self\):
    &quot;&quot;&quot; classify the data into land, coast, missing, hole &quot;&quot;&quot;
    data = self\.data\[&#39;data&#39;\]
    self\.header = self\.data\[&#39;header&#39;\]

    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&gt;250\] = 0
    self\.ice = data

def geocoordinate\(self\):
    &quot;&quot;&quot; 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\[&#39;ROWS&#39;\]
        COLS = self\.header\[&#39;COLS&#39;\]
    except:
        raise AttributeError\(&#39;object needs to have header, \\
                did you run self\.classify?&#39;\)

    datadir = &#39;nsidc0081\_ssmi\_nrt\_seaice&#39;

    lonfile = os\.path\.join\(datadir,&#39;psn25lons\_v2\.dat&#39;\)
    lons = np\.fromfile\(lonfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/100000\.
    self\.lons = lons\.reshape\(COLS,ROWS\)

    latfile = os\.path\.join\(datadir,&#39;psn25lats\_v2\.dat&#39;\)
    lats = np\.fromfile\(latfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/100000\.
    self\.lats = lats\.reshape\(COLS,ROWS\)

    areafile = os\.path\.join\(datadir,&#39;psn25area\_v2\.dat&#39;\)
    area = np\.fromfile\(latfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/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

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

···

On 10/25/10 2:13 PM, John wrote:

http://picasaweb.google.com/lh/photo/6Dnylo0BcjX0A-wdmczmNg?feat=directlink

Any ideas??

-john

On Mon, Oct 25, 2010 at 9:53 PM, John<washakie@...287...> 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<washakie@...287...> 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<jswhit@...146...> 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.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

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

···

On Mon, Oct 25, 2010 at 10:24 PM, Jeff Whitaker <jswhit@...146...> 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<washakie@...287...> 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<washakie@...287...> 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=&#39;stere&#39;,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<jswhit@...146...> >>>> 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\):
    &quot;&quot;&quot; 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
    &quot;&quot;&quot;
    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 \{&#39;header&#39;:header,&#39;data&#39;:seaiceconc\}

def conv2percentage\(self\):
    self\.seaicepercentage = self\.data\[&#39;data&#39;\]/2\.5

def classify\(self\):
    &quot;&quot;&quot; classify the data into land, coast, missing, hole &quot;&quot;&quot;
    data = self\.data\[&#39;data&#39;\]
    self\.header = self\.data\[&#39;header&#39;\]

    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&gt;250\] = 0
    self\.ice = data

def geocoordinate\(self\):
    &quot;&quot;&quot; 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\[&#39;ROWS&#39;\]
        COLS = self\.header\[&#39;COLS&#39;\]
    except:
        raise AttributeError\(&#39;object needs to have header, \\
                did you run self\.classify?&#39;\)

    datadir = &#39;nsidc0081\_ssmi\_nrt\_seaice&#39;

    lonfile = os\.path\.join\(datadir,&#39;psn25lons\_v2\.dat&#39;\)
    lons = np\.fromfile\(lonfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/100000\.
    self\.lons = lons\.reshape\(COLS,ROWS\)

    latfile = os\.path\.join\(datadir,&#39;psn25lats\_v2\.dat&#39;\)
    lats = np\.fromfile\(latfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/100000\.
    self\.lats = lats\.reshape\(COLS,ROWS\)

    areafile = os\.path\.join\(datadir,&#39;psn25area\_v2\.dat&#39;\)
    area = np\.fromfile\(latfile,dtype=np\.dtype\(&#39;i4&#39;\)\)/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.Whitaker@...2583...9...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

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

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 10/25/10 2:28 PM, John wrote:

On Mon, Oct 25, 2010 at 10:24 PM, Jeff Whitaker<jswhit@...146...> 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<washakie@...287...> 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<washakie@...287...> 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<jswhit@...146...> >>>>> 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.Whitaker@...259...
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.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Jeff and others,

I have a working solution and a non-working solution. I'm wondering
whether someone can help me understand the problems in the latter.
Please forgive in advance the long post, but all of my code is here.

First, you can get a dataset from :
ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/

and the projection files here:
http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats

I've made a class for handling this data:
<code>
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 pflexpart 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']))
        seaiceconc = seaiceconc.T

        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 =
'/flex_wrk/jfb/RESEARCH_ARCTIC/SEAICE/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(ROWS,COLS)

        latfile = os.path.join(datadir,'psn25lats_v2.dat')
        lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
        self.lats = lats.reshape(ROWS,COLS)

        areafile = os.path.join(datadir,'psn25area_v2.dat')
        area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
        self.area = area.reshape(ROWS,COLS)

    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=-45,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[1]
        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,0.5)
        lats = np.arange(-90,90,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, xin, yin,\
                  xout, yout, masked=True)

        x,y = np.meshgrid(xin,yin) # xin, yin are 1-d, x,y are 2-d
(same shape as data)
        map.contourf(x,y,np.flipud(self.ice),clevs=np.arange(0,250,1))

        self.iceglobe = dataout
        self.m = map

</code>

You can use this as follows:

import NSIDC
d = NSIDC('nt_20080827_f17_nrt_n.bin')
d.classify()
d.geocoordinate()
d.regrid2globe()
d.m.drawcoastlines()

This should produce a plot.

QUESTION ONE: It seems possibly I'm slightly off (look at the great
lakes). Any suggestions as to why?

QUESTION TWO: Please, suggested improvements, code review or
simplification is welcome.

Now, on to the NON-WORKING example. I have a function (it's not too
pretty, but it basically takes a bunch of pre-defined regions and
plots a 'grid' of data to them. I've tried to make the necessary edits
so you can use the code as is.

data = d.iceglobe # using the d from the above example
lons = np.arange(-180,180,0.5)
lats = np.arange(-90,90,0.5)

fig,m = plot_grid((lons,lats,data),region=None)

QUESTION THREE: Notice that I'm quite sure the grid is upside down
(you can see again the great lakes and Hudson bay outline in the top
left). However, I've tried a variety of np.fliplr / np.flipud / data.T
combinations, but I can't seem to get it projected properly. Any ideas
here??

QUESTION FOUR: As above in Q2.

<code>

def plot_grid(D,region=None,dres=0.5,
              transform=True,figname=None,fillcontinents=False,
              points=False):
    """ plot an array over a basemap. The required argument "D" is
    either:
    * a tuple of (x,y,z)
    * a tuple with (lon,lat,grid)
    * a tuple with (grid,)

    Usage::

        > FIG = plot_grid(D,**kwargs)

    =============== ========================================
    keyword description
    =============== ========================================
    dres resolution of the grid
    fillcontinents fills continents if True
    plotargs dictionary of plot arguments
    figname A figurename, if passed will save the
                        figure with figname
    points set to True if passing a x,y,z matrix of
                        points
    region A region from :func:`get_base1`
    =============== ========================================

    """
    print "length of D: %s" % len(D)

    if isinstance(D,np.ndarray):
        assert len(D.shape) == 2, "D grid must be 2d"
    #if len(D)==1:
        #assume full earth grid with dres
        print 'received grid of shape:', D.shape
        if D.shape[0] == 720:
            lons = np.arange(-180,180,dres)
        elif D.shape[0] == 721:
            lons = np.arange(-180,180.01,dres)
        if D.shape[1] == 360:
            lats = np.arange(-90,90,dres)
        elif D.shape[1] == 361:
            lats = np.arange(-90,90.01,dres)
        points = False
        z = D.T

    elif len(D)==3:
        points = True
        x = D[0]
        y = D[1]
        z = D[2]

        if len(z.shape) > 1:
            points = False
            lons = x
            lats = y

    #Set up a basemap

    ## CHANGED THIS HERE FOR THE CODEX ##
    ## Be sure to set region=None ##
    if isinstance(region,str):
        print "getting basemap with region: %s" % region
        fig,m = get_base1(region=region)
    else:
        fig = plt.figure()
        ax = fig.add_axes()
        a = 6378.273e3
        ec = 0.081816153
        b = a*np.sqrt(1.-ec**2)
        m = Basemap(projection='stere',lat_0=90,lon_0=-45,lat_ts=70,\
                      llcrnrlat=33.92,llcrnrlon=279.96,\
                      urcrnrlon=102.34,urcrnrlat=31.37,\
                      rsphere=(a,b))

        m.drawcoastlines()

    if points == True:
    # Plot individual data points
        norm = colors.normalize(z.min(),z.max())
        for i in range(len(y)):
            xpt,ypt = m(x[i],y[i])
            cmap = cm.jet(norm(z[i]))
            #cmap = 'k'
            m.plot([xpt],[ypt],'.',color=cmap,markersize=2)
            #plt.plot(x[i],y[i],'.',color=cmap,markersize=2)

    if points == False:
        #transform Z data into projection
        #transform to nx x ny regularly spaced native projection grid
        dx = 2.*np.pi*m.rmajor/len(lons)
        nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1
        if transform:
            # Need this if we choose lon,lat approach
            Zt,xx,yy = m.transform_scalar(z,lons,lats,nx,ny,returnxy=True)
        else:
            Zt = z
        print m.projection

        m.imshow(Zt)
        if fillcontinents:
            m.fillcontinents()
        plt.colorbar()
        #plt.imshow(Z)
    if figname != None:
        #plt.ylim([40,90])
        #plt.title('data locations on mercator grid')
        plt.savefig(figname)
    else:
        plt.show()

</code>