Average seaice 1994 to 2000 during Spring

Create contourmaps of sea-ice using basemap

Average seaice 1994 to 2000 during Spring
Average seaice 1994 to 2000 during Spring

Often I need to visualize non-uniform spatially distributed geographic data. I have found that a combination of Numpy and Basemap creates a powerful and efficient combination of tools that produce print quality maps for visualizing my data.  I ususally choose to interpolate my data to a grid using Numpy which I then plot using filled contours in Basemap.  The image to the left shows an example of how I use Numpy and Basemap to create a map of the average spring sea-ice extent between 1994 and 2000. The data used here are from the Walsh and Chapman sea-ice database, which is uniformly spaced at 1 x 1 degrees latitude/longitude. If you wanted to use non-uniform data then you just have to interpolate your data to a grid prior to plotting.  The interpolation is done using interp2 and the grid is creating using meshgrid.

To get started making a map of sea-ice download the data from here. When I started using these data I wanted to show the seasonal variability in sea-ice extent between years, so I wrote a Python script that estimates the seasonal average sea-ice concentration mean_ice.py (862).  The output of mean_ice.py (862) is then used to create a contourfilled map of the ice concentration over the North Pole.  Note that these data do not cover the South Pole, but there are extensive satellite databases that do.  The map is created using the script contourICEMaps.py (950).

Note that these scripts require a number of Python modules including the Python NetCDF4 interface, Basemap, Numpy, and Matplotlib. Details on how to install all these modules are found here.

download mean_ice.py (862)

import os, sys, math, string
import numpy as np
# Program to calculate the average of ice concentration over seasons (months).
# This proggram calculates the average ice concentration for the whole domain
# for a given season and range of years. For example, for each season in the
# list seasons is looped over and concentration extracted for the period defined by
# setting startYear and endYear. The result is written to file and the result can
# be plotted using contourMaps.py.
#
# Trond Kristiansen, 03.10.03, revised 08.03.2011, 09.03.2011 at Kens house

seasons=["winter","spring","summer","autumn"]

startYear = 1994
endYear = 2000 # note: endyear is not included

years=[year for year in range(startYear,endYear)]

for season in seasons:

print "\n Extracting ice data for season %s"%(season)
# Check for which season to average over
if season == 'winter':
months = [1,2,3]
if season == 'spring':
months = [4,5,6]
if season == 'summer':
months = [7,8,9]
if season == 'autumn':
months = [10,11,12]

# Read the first header
file = open('arctic.latlon','r')

arrayLat=[]
arrayLon=[]

# read the first 58 lines, which is the latitude matrix into
# an array. The matrix is 80X58 and covers the entire northern
# hemisphere.
for i in range(58):
list=[]
line = file.readline()
if not line:
break
l=string.split(line)
for j in range(len(l)):
list.append(l[j])

arrayLat.append(list)

# Read the 58 last lines which is the longitude matrix into
# an array.
for i in range(58):
list2=[]
line2 = file.readline()
if not line2:
break
l2=string.split(line2)
for j in range(len(l2)):
list2.append(l2[j])

arrayLon.append(list2)

#print '*** TEST***'
#print 'Compare to the README file of JOHN WALSH for testing the output'

#print '(1,1) = %sN %sE'%(arrayLat[0][0],arrayLon[0][0])
#print '(1,2) = %sN %sE'%(arrayLat[1][0],arrayLon[1][0])
#print '(35,24)= %sN %sE'%(arrayLat[23][34],arrayLon[23][34])
#print '(1,58) = %sN %sE'%(arrayLat[57][0],arrayLon[57][0])
#print '(80,1) = %sN %sE'%(arrayLat[0][79],arrayLon[0][79])
#print '(80,58)= %sN %sE'%(arrayLat[57][79],arrayLon[57][79])

# The result should be
# (1,1) = 048.9512N 125.9228E
# (1,2) = 049.5031N 127.0948E
# (35,24)= 090.0000N 250.0000E
# (1,58) = 041.9167N 205.0000E
# (80,1) = 039.4629N 007.0721E
# (80,58)= 033.5996N 302.9269E

file.close()

# Continue by reading the data.
# Now we have the array contining the lat and long of
# each grid point.

# starting with the first header

icematrix=np.zeros((len(years),58,80))
y=0

for year in years:

file = open('icecon.01-99','r')
line = file.readline()

while line:
l = string.split(line)
month_read = l[0]
year_read = l[1]

if (int(year_read) == int(year)) and (int(month_read) in months):
print 'Reading data for month %s and year %s'%(month_read,year_read)
j=0
latitude=[]
longitude=[]
lineNumber=0
while int(j) < 58:
line = file.readline()

if not line:
break
i=0
while int(i) < 80: ice = line[i] #print line #print line[i], i latitude = arrayLat[lineNumber][i] longitude = arrayLon[lineNumber][i] if ice == '*': ice=10 if ice == '.': ice='land' if float(longitude) > 180:
longitude= float(longitude) - 360
if (str(ice) !='land'):
icematrix[y,j,i]=icematrix[y,j,i] + int(ice)
#print "icematrix",icematrix[y,j,i]
i+=1
lineNumber=lineNumber+1
j+=1

# End of the month, read the next header (month-year)
line = file.readline()
if not line:break
l = string.split(line)
month_read = l[0]
year_read = l[1]
else:
# Keep reading month for month until a selected one occurs
# one month = 58*81 bytes
file.seek(58*81,1)
line = file.readline()
if not line:break
l = string.split(line)
month_read = l[0]
year_read = l[1]
#print "Checking ",month_read, year_read
y+=1

file.close()

# Now calculate the average ice concentration for the season at each grid point.

FileOut="ICE_average_"+str(startYear)+"_to_"+str(endYear)+'_'+str(season)+'.txt'
print 'Output file is named: %s'%(FileOut)
if os.path.exists(FileOut):
os.remove(FileOut)
f = open(FileOut,'a')

print 'Averaging ice-concentration in season %s in year %s\n'%(season,year)

for j in range(58):
for i in range(80):
latitude = arrayLat[j][i]
longitude = arrayLon[j][i]

mean_ice=np.mean((icematrix[:,j,i]/float(len(months))))
data = '%.3f\t%.3f\t%.3f\n'%(float(longitude),float(latitude),mean_ice)

f.writelines(data)
i=0

f.close()

download contourICEMaps.py (950) Laplace Filter (1160) mpl_util.py (1091)

import os, sys, datetime, string
import numpy as np
from netCDF4 import Dataset
import numpy.ma as ma
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, shiftgrid, NetCDFFile
from pylab import *
import laplaceFilter
import mpl_util

__author__ = 'Trond Kristiansen'
__email__ = 'trond.kristiansen (at) imr.no'
__created__ = datetime.datetime(2008, 8, 15)
__modified__ = datetime.datetime(2009, 7, 21)
__version__ = "1.0"
__status__ = "Development"

def findSubsetIndices(min_lat,max_lat,min_lon,max_lon,lats,lons):
"""Array to store the results returned from the function"""
res=np.zeros((4),dtype=np.float64)
minLon=min_lon; maxLon=max_lon
distances1 = []; distances2 = []
indices=[]; index=1
for point in lats:
s1 = max_lat-point # (vector subtract)
s2 = min_lat-point # (vector subtract)
distances1.append((np.dot(s1, s1), point, index))
distances2.append((np.dot(s2, s2), point, index-1))
index=index+1
distances1.sort()
distances2.sort()
indices.append(distances1[0])
indices.append(distances2[0])
distances1 = []; distances2 = []; index=1

for point in lons:
s1 = maxLon-point # (vector subtract)
s2 = minLon-point # (vector subtract)
distances1.append((np.dot(s1, s1), point, index))
distances2.append((np.dot(s2, s2), point, index-1))
index=index+1
distances1.sort()
distances2.sort()
indices.append(distances1[0])
indices.append(distances2[0])

""" Save final product: max_lat_indices,min_lat_indices,max_lon_indices,min_lon_indices"""
minJ=indices[1][2]
maxJ=indices[0][2]
minI=indices[3][2]
maxI=indices[2][2]
res[0]=minI; res[1]=maxI; res[2]=minJ; res[3]=maxJ;
return res

def makeMap(lonStart,lonEnd,latStart,latEnd,name,stLon,stLat):
plt.figure(figsize=(8,8))
"""Get the etopo2 data"""
etopo1name='ETOPO1_Ice_g_gmt4.grd'
etopo1 = Dataset(etopo1name,'r')
lons = etopo1.variables["lon"][:]
lats = etopo1.variables["lat"][:]

res = findSubsetIndices(latStart-5,latEnd+5,lonStart-40,lonEnd+10,lats,lons)
lon,lat=np.meshgrid(lons[res[0]:res[1]],lats[res[2]:res[3]])
print "Extracted data for area %s : (%s,%s) to (%s,%s)"%(name,lon.min(),lat.min(),lon.max(),lat.max())
bathy = etopo1.variables["z"][int(res[2]):int(res[3]),int(res[0]):int(res[1])]
bathySmoothed = laplaceFilter.laplace_filter(bathy,M=None)
levels=[-6000,-5000,-3000, -2000, -1500, -1000,-500, -400, -300, -250, -200, -150, -100, -75, -65, -50, -35, -25, -15, -10, -5, 0]
if lonStart< 0 and lonEnd < 0:
lon_0= - (abs(lonEnd)+abs(lonStart))/2.0
else:
lon_0=(abs(lonEnd)+abs(lonStart))/2.0
map = Basemap(llcrnrlat=latStart,urcrnrlat=latEnd,\
llcrnrlon=lonStart,urcrnrlon=lonEnd,\
rsphere=(6378137.00,6356752.3142),\
resolution='l',area_thresh=1000.,projection='lcc',\
lat_1=latStart,lon_0=lon_0)
x, y = map(lon,lat)
map.drawcoastlines()
map.drawcountries()
map.fillcontinents(color='grey')
map.drawmeridians(np.arange(lons.min(),lons.max(),10),labels=[0,0,0,1])
map.drawparallels(np.arange(lats.min(),lats.max(),4),labels=[1,0,0,0])
#map.bluemarble()

CS1 = map.contourf(x,y,bathySmoothed,levels,
cmap=mpl_util.LevelColormap(levels,cmap=cm.Blues_r),
extend='upper',
alpha=1.0,
origin='lower')
CS1.axis='tight'
"""Plot the station as a position dot on the map"""
pt,ypt = map(stLon,stLat)
map.plot([xpt],[ypt],'ro', markersize=10)
plt.text(xpt+100000,ypt+100000,name)
plt.title('Area %s'%(name))
plotfile='figures/map_'+str(name)+'.pdf'
plt.savefig(plotfile,dpi=150,orientation='portrait')
plt.show()

def main():
names=['bayOfBiscay']
lat_start=[40]
lat_end =[81]
lon_start=[-10]
lon_end =[75]
"""List of stations for each area"""
stationlonlist=[ 2.4301, -22.6001, -47.0801, 13.3801, -67.2001]
stationlatlist=[54.5601, 63.7010, 60.4201, 67.5001, 41.6423]
for i in range(len(lat_start)):
if i==0:
makeMap(lon_start[i],lon_end[i],lat_start[i],lat_end[i],names[i],stationlonlist[i],stationlatlist[i])
if __name__ == "__main__":
main()

Leave a Reply

Your email address will not be published. Required fields are marked *