name = '2015-04-20-arctic_sea_ice_concentration'
title = "Arctic Sea Ice Concentration"
import os
from datetime import datetime
from IPython.core.display import HTML
with open('creative_commons.txt', 'r') as f:
html = f.read()
html = '''
<small>
<p> This post was written as an IPython notebook.
It is available for <a href='https://ocefpaf.github.com/python4oceanographers/downloads/notebooks/%s.ipynb'>download</a>
or as a static <a href='https://nbviewer.ipython.org/url/ocefpaf.github.com/python4oceanographers/downloads/notebooks/%s.ipynb'>html</a>.</p>
<p></p>
%s''' % (name, name, html)
%matplotlib inline
from matplotlib import style
style.use('ggplot')
hour = datetime.utcnow().strftime('%H:%M')
comments="true"
date = '-'.join(name.split('-')[:3])
slug = '-'.join(name.split('-')[3:])
metadata = dict(title=title,
date=date,
hour=hour,
comments=comments,
slug=slug,
name=name)
markdown = """Title: {title}
date: {date} {hour}
comments: {comments}
slug: {slug}
{{% notebook {name}.ipynb cells[1:] %}}
""".format(**metadata)
content = os.path.abspath(os.path.join(os.getcwd(), os.pardir, os.pardir, '{}.md'.format(name)))
with open('{}'.format(content), 'w') as f:
f.writelines(markdown)
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.
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.)
!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:
dx = dy = 25000
x = np.arange(-3850000, +3750000, +dx)
y = np.arange(+5850000, -5350000, -dy)
x.shape, y.shape, ice.shape
((304,), (448,), (448, 304))
Finally the figure!
%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))
HTML(html)
This post was written as an IPython notebook. It is available for download or as a static html.