|
Size: 982
Comment:
|
← Revision 7 as of 2008-12-10 12:54:38 ⇥
Size: 5993
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 12: | Line 12: |
| * Download netcdf ice concentration data [[ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/psi-concentration/data|ftp server]] in [[http://nsidc.org/data/grids/ps_grid.html|polar stereographic projection]] * Read in the data using [[http://www.pyngl.ucar.edu/Nio.shtml|PyNIO]] |
* Download ice concentration data [[ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/psi-concentration/data|ftp server]] * Read in the netcdf data using [[http://www.pyngl.ucar.edu/Nio.shtml|PyNIO]] |
| Line 15: | Line 15: |
| * 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 [[http://nsidc.org/data/grids/ps_grid.html|polar stereographic projection]] |
* 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 [[http://nsidc.org/data/grids/ps_grid.html|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 {{{#!python from numpy import * # Polar stereographic (NSDIDC grid) def mapxy(x, y, sgn): cdr = 57.29577951 slat = 70.0 re = 6378.273 ec2 = .006693883 ec = sqrt(ec2) if sgn == 1: delta = 45.0 else: delta = 0.0 sl = slat * pi/180.0 rho = sqrt(x**2 + y**2) cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2)) t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0) if (absolute(slat-90.0) < 1.e-5): t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re else: t = rho * t / (re * cm) chi = (pi / 2.0) - 2.0 * arctan(t) 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) alat = sgn * alat along = arctan2(sgn * x, -sgn * y) along = sgn * along along = along * 180. / pi alat = alat * 180. / pi along = along - delta return [alat,along] def mapll(lat, lon,sgn): cdr = 57.29577951 slat = 70. re = 6378.273 ec2 = .006693883 ec = sqrt(ec2) if sgn == 1: delta = 45.0 else: delta = 0.0 latitude = absolute(lat) * pi/180. longitude = (lon + delta) * pi/180. t = tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2) if ((90-slat) < 1.e-5): rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec)) else: sl = slat*pi/180. tc = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.) mc = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2)) rho = re*mc*t/tc y=-rho*sgn*cos(sgn*longitude) x= rho*sgn*sin(sgn*longitude) return [x,y] }}} === Solution === {{{#!python # imports from numpy import * from scipy import array,arange,nan from PyNGL import Ngl from PyNGL import Nio # Polar stereographic (NSDIDC grid) def mapxy(x, y, sgn): cdr = 57.29577951 slat = 70.0 re = 6378.273 ec2 = .006693883 ec = sqrt(ec2) if sgn == 1: delta = 45.0 else: delta = 0.0 sl = slat * pi/180.0 rho = sqrt(x**2 + y**2) cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2)) t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0) if (absolute(slat-90.0) < 1.e-5): t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re else: t = rho * t / (re * cm) chi = (pi / 2.0) - 2.0 * arctan(t) 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) alat = sgn * alat along = arctan2(sgn * x, -sgn * y) along = sgn * along along = along * 180. / pi alat = alat * 180. / pi along = along - delta return [alat,along] def mapll(lat, lon,sgn): cdr = 57.29577951 slat = 70. re = 6378.273 ec2 = .006693883 ec = sqrt(ec2) if sgn == 1: delta = 45.0 else: delta = 0.0 latitude = absolute(lat) * pi/180. longitude = (lon + delta) * pi/180. t = tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2) if ((90-slat) < 1.e-5): rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec)) else: sl = slat*pi/180. tc = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.) mc = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2)) rho = re*mc*t/tc y=-rho*sgn*cos(sgn*longitude) x= rho*sgn*sin(sgn*longitude) return [x,y] def mapij(x,y,sgn): #i=(x+3850)/12.5 #j=(abs(y)+5850)/12.5 #ir=round((x+3850)/12.5) #jr=round((abs(y)+5850)/12.5) i=(x+5850)/12.5 j=(abs(y)+3850)/12.5 ir=round((x+5850)/12.5) jr=round((abs(y)+3850)/12.5) return [i,j,ir,jr] # Main #nc = Nio.open_file('20080915.nc') nc = Nio.open_file('20081122.nc') C=array(nc.variables['concentration'])[0,:,:].astype(float) C1=array(nc.variables['concentration'])[0,:,:].astype(float) C1.fill(128.0) latList = linspace(75.0, 85.0, 100) resultList=[] print nc for i in latList: stereoCoords=mapll(i,0.0,1) print stereoCoords gridCoords=mapij(stereoCoords[0],stereoCoords[1],1) print gridCoords if C[gridCoords[2],gridCoords[3]] != 0.0 and C[gridCoords[2],gridCoords[3]] != 128.0 : resultList.append([i, C[gridCoords[2],gridCoords[3]],gridCoords[2],gridCoords[3]]) C1[gridCoords[2],gridCoords[3]]=20.0 resultList.sort() print resultList[0] }}} |
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]
Solution
1 # imports
2 from numpy import *
3 from scipy import array,arange,nan
4 from PyNGL import Ngl
5 from PyNGL import Nio
6
7 # Polar stereographic (NSDIDC grid)
8 def mapxy(x, y, sgn):
9 cdr = 57.29577951
10 slat = 70.0
11 re = 6378.273
12 ec2 = .006693883
13 ec = sqrt(ec2)
14 if sgn == 1:
15 delta = 45.0
16 else:
17 delta = 0.0
18 sl = slat * pi/180.0
19 rho = sqrt(x**2 + y**2)
20 cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2))
21 t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0)
22 if (absolute(slat-90.0) < 1.e-5):
23 t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re
24 else:
25 t = rho * t / (re * cm)
26 chi = (pi / 2.0) - 2.0 * arctan(t)
27 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 *
28 chi)+ (7.0 * ec2**3.0 / 120.0) * sin(6.0 * chi)
29 alat = sgn * alat
30 along = arctan2(sgn * x, -sgn * y)
31 along = sgn * along
32
33 along = along * 180. / pi
34 alat = alat * 180. / pi
35 along = along - delta
36 return [alat,along]
37
38 def mapll(lat, lon,sgn):
39 cdr = 57.29577951
40 slat = 70.
41 re = 6378.273
42 ec2 = .006693883
43 ec = sqrt(ec2)
44 if sgn == 1:
45 delta = 45.0
46 else:
47 delta = 0.0
48 latitude = absolute(lat) * pi/180.
49 longitude = (lon + delta) * pi/180.
50 t = tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2)
51 if ((90-slat) < 1.e-5):
52 rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec))
53 else:
54 sl = slat*pi/180.
55 tc = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.)
56 mc = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2))
57 rho = re*mc*t/tc
58
59 y=-rho*sgn*cos(sgn*longitude)
60 x= rho*sgn*sin(sgn*longitude)
61 return [x,y]
62
63 def mapij(x,y,sgn):
64 #i=(x+3850)/12.5
65 #j=(abs(y)+5850)/12.5
66 #ir=round((x+3850)/12.5)
67 #jr=round((abs(y)+5850)/12.5)
68 i=(x+5850)/12.5
69 j=(abs(y)+3850)/12.5
70 ir=round((x+5850)/12.5)
71 jr=round((abs(y)+3850)/12.5)
72 return [i,j,ir,jr]
73
74 # Main
75 #nc = Nio.open_file('20080915.nc')
76 nc = Nio.open_file('20081122.nc')
77 C=array(nc.variables['concentration'])[0,:,:].astype(float)
78 C1=array(nc.variables['concentration'])[0,:,:].astype(float)
79 C1.fill(128.0)
80 latList = linspace(75.0, 85.0, 100)
81 resultList=[]
82 print nc
83 for i in latList:
84 stereoCoords=mapll(i,0.0,1)
85 print stereoCoords
86 gridCoords=mapij(stereoCoords[0],stereoCoords[1],1)
87 print gridCoords
88 if C[gridCoords[2],gridCoords[3]] != 0.0 and C[gridCoords[2],gridCoords[3]] != 128.0 :
89 resultList.append([i, C[gridCoords[2],gridCoords[3]],gridCoords[2],gridCoords[3]])
90 C1[gridCoords[2],gridCoords[3]]=20.0
91 resultList.sort()
92 print resultList[0]