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 ice concentration data ftp server
Read in the netcdf data using PyNIO
To calculate the meridional ice concentration do the following:
- Loop over latitude 75°N to 85°N
- 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

X,Y=mapll(lat,lon,sgn) converts from latitude and longitude to X,Y
lat,lon=mapxy(X,Y,sgn) converts from X,Y [km] to latitude and longitude
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]