Recently I got a request to help with this question. I decided to create a post out of it because it is a nice example to show how to do something from start to finish. (Something that I have to do for my SWC training, but that is another story.)

First download the data from http://nsidc.org/data/NSIDC-0081. The binary file description says:

> One-byte flat scaled binary (preceded by a 300-byte header)

that means we can read the data with numpy's fromfile after skipping the header. We also need to re-size the matrix to the grid specifications:

> North: 304 columns x 448 rows

According to the documentation the sea ice concentration is scaled by 250, to make it simpler to plot let's just scaled it to 0-1. The final touch is to mask out the land values.

In [2]:

import numpy as np
import numpy.ma as ma


infile = './data/nt_20150326_f17_nrt_n.bin'
with open(infile, 'rb') as fr:
    hdr = fr.read(300)
    ice = np.fromfile(fr, dtype=np.uint8)

ice = ice.reshape(448, 304)

ice = ice / 250.

ice = ma.masked_greater(ice, 1.0)

The tricky part is to find the right projection to plot the data. Luckily, the text file that accompanies the data has everything we need. Including a proj4 string! (In case you are using Basemap or a future version of cartopy that is all you will need.)

In [3]:

!cat ./data/nt_20150326_f17_nrt_n.bin.txt

Name: nt_20150326_f17_nrt_n.bin
Grid Name: Polar Stereo Northern 25km
NSIDC GPD File Name: ftp://sidads.colorado.edu/pub/tools/mapx/nsidc_maps/N3B.gpd
Data Offset: 300
Data Type: UNSIGNED_INTEGER
Bytes Per Pixel:1
Byte Order: n/a
Image Width: 304
Image Height: 448
Proj.4 Projection: +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs
Missing Data Value: 255
Upper Left Corner X Coordinate: -3850000.0
Upper Left Corner Y Coordinate: 5850000.0
Lower Right Corner X Coordinate: 3750000.0
Lower Right Corner Y Coordinate: -5350000.0
Subset Parent Grid Cell Start Row: 0.0
Subset Parent Grid Cell End Row: 448.0
Subset Parent Grid Cell Start Column: 0.0
Subset Parent Grid Cell End Column: 304.0
Minimum Value: 0.0
Maximum Value: 255.0
History:

In [4]:

dx = dy = 25000

x = np.arange(-3850000, +3750000, +dx)
y = np.arange(+5850000, -5350000, -dy)

x.shape, y.shape, ice.shape

Out[4]:

((304,), (448,), (448, 304))

Finally the figure!

In [5]:

%matplotlib inline

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(9, 9))
ax = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0))

cs = ax.coastlines(resolution='110m', linewidth=0.5)

ax.gridlines()
ax.set_extent([-180, 180, 40, 90], crs=ccrs.PlateCarree())

kw = dict(central_latitude=90, central_longitude=-45, true_scale_latitude=70)
cs = ax.pcolormesh(x, y, ice, cmap=plt.cm.Blues,
                   transform=ccrs.Stereographic(**kw))

In [6]:

HTML(html)

Out[6]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

Arctic Sea Ice Concentration

Comments