Data formats

Netcdf

What is the northernmost ice free position on a given longitude? Calculate the latitude for the prime meridian on the 15th September 2008.

Download netcdf ice concentration data ftp server in polar stereographic projection
Read in the data using PyNIO
To calculate the meridional ice concentration do the following:
Loop over latitude 75N to 85N
Convert to polar stereographic coordinates [km].
Convert to pixel indices (y,x). The resolution of the grid is 12.5 km. The offset to the pole is given in the definition polar stereographic projection

Polar stereographic projection

lat,lon=mapll(X,Y,sgn) converts from X,Y [km] to latitude and longitude
X,Y=mapxy(lat,lon,sgn) converts from latitude and longitude to X,Y
Set sgn=1 for north and sgn=-1 for southern hemisphere

   1 from numpy import *
   2 
   3 # Polar stereographic (NSDIDC grid)
   4 def mapxy(x, y, sgn):
   5     cdr  = 57.29577951
   6     slat = 70.0
   7     re   = 6378.273
   8     ec2   = .006693883
   9     ec    =  sqrt(ec2)
  10     if sgn == 1:
  11         delta = 45.0
  12     else:
  13         delta = 0.0
  14     sl = slat * pi/180.0
  15     rho = sqrt(x**2 + y**2)
  16     cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2))
  17     t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0)
  18     if  (absolute(slat-90.0) < 1.e-5):
  19         t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re
  20     else:
  21         t = rho * t / (re * cm)
  22     chi = (pi / 2.0) - 2.0 * arctan(t)
  23     alat = chi + ((ec2 / 2.0) + (5.0 * ec2**2.0 / 24.0) + (ec2**3.0 / 12.0)) * sin(2 * chi) + ((7.0 * ec2**2.0 / 48.0) + (29.0 * ec2**3 / 240.0)) *sin(4.0 * chi)+ (7.0 * ec2**3.0 / 120.0) * sin(6.0 * chi)
  24     alat = sgn * alat
  25     along = arctan2(sgn * x, -sgn * y)
  26     along = sgn * along
  27 
  28     along = along * 180. / pi
  29     alat  = alat * 180. / pi
  30     along = along - delta
  31     return [alat,along]
  32 
  33 def  mapll(lat, lon,sgn):
  34     cdr  = 57.29577951
  35     slat = 70.
  36     re   = 6378.273
  37     ec2   = .006693883
  38     ec    =  sqrt(ec2)
  39     if sgn == 1:
  40         delta = 45.0
  41     else:
  42         delta = 0.0
  43     latitude  = absolute(lat) * pi/180.
  44     longitude = (lon + delta) * pi/180.
  45     t =   tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2)
  46     if ((90-slat) < 1.e-5):
  47         rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec))
  48     else:
  49         sl  = slat*pi/180.
  50         tc  = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.)
  51         mc  = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2))
  52         rho = re*mc*t/tc
  53 
  54     y=-rho*sgn*cos(sgn*longitude)
  55     x= rho*sgn*sin(sgn*longitude)
  56     return [x,y]