|
Size: 3053
Comment:
|
Size: 3053
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 25: | Line 25: |
| * {{{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}}} |
* {{{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 |
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]