I have ported from Fortran to Python an algorithm that calculates the numerical derivative along the x direction (longitudinal) of a scalar function s on a rectilinear grid that has equal grid spacing in the x and y directiondirections (2.5 degrees x 2.5 degrees). The input data s is oriented south to north for the latitude and west to east for the longitude. The order of the derivatives should be in the following way:
At the poles, the derivatives isare set to zero, i.e., dsdx[j,0] = 0 and dsdx[j,-1] = 0.
Both lat and lon are one dimensional-dimensional arrays, having such values:
[-90. -87.5 -85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5 -35. -32.5
-30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5 -10. -7.5 -5. -2.5
0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5 25. 27.5
30. 32.5 35. 37.5 40. 42.5 45. 47.5 50. 52.5 55. 57.5
60. 62.5 65. 67.5 70. 72.5 75. 77.5 80. 82.5 85. 87.5
90. ]
[-90. -87.5 -85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5 -35. -32.5
-30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5 -10. -7.5 -5. -2.5
0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5 25. 27.5
30. 32.5 35. 37.5 40. 42.5 45. 47.5 50. 52.5 55. 57.5
60. 62.5 65. 67.5 70. 72.5 75. 77.5 80. 82.5 85. 87.5
90. ]
[ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5
25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5
50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5
75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5
100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5
125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5
150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5
175. 177.5 -180. -177.5 -175. -172.5 -170. -167.5 -165. -162.5
-160. -157.5 -155. -152.5 -150. -147.5 -145. -142.5 -140. -137.5
-135. -132.5 -130. -127.5 -125. -122.5 -120. -117.5 -115. -112.5
-110. -107.5 -105. -102.5 -100. -97.5 -95. -92.5 -90. -87.5
-85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5
-35. -32.5 -30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5
-10. -7.5 -5. -2.5]
[ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5
25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5
50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5
75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5
100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5
125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5
150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5
175. 177.5 -180. -177.5 -175. -172.5 -170. -167.5 -165. -162.5
-160. -157.5 -155. -152.5 -150. -147.5 -145. -142.5 -140. -137.5
-135. -132.5 -130. -127.5 -125. -122.5 -120. -117.5 -115. -112.5
-110. -107.5 -105. -102.5 -100. -97.5 -95. -92.5 -90. -87.5
-85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5
-35. -32.5 -30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5
-10. -7.5 -5. -2.5]
8.50002
8.8
9.00002
9.20001
9.50002
9.70001
9.8
10.0
10.1
10.2
10.3
10.4
10.5
10.5
10.5
10.5
10.5
10.5
10.4
10.3
10.2
10.1
10.0
9.8
9.60001
9.40001
9.20001
9.00002
8.70001
8.50002
8.20001
7.90001
7.60001
7.3
6.90001
6.60001
8.50002
8.8
9.00002
9.20001
9.50002
9.70001
9.8
10.0
10.1
10.2
10.3
10.4
10.5
10.5
10.5
10.5
10.5
10.5
10.4
10.3
10.2
10.1
10.0
9.8
9.60001
9.40001
9.20001
9.00002
8.70001
8.50002
8.20001
7.90001
7.60001
7.3
6.90001
6.60001
8.99290394761e-07
7.19448782308e-07
5.39579725689e-07
3.59710669071e-07
1.79869056618e-07
0.0
-1.79869056618e-07
-1.79869056618e-07
-5.39579725689e-07
-7.19421338143e-07
-7.19421338143e-07
-7.19421338143e-07
-8.99290394761e-07
-1.07915945138e-06
-1.07915945138e-06
-1.07915945138e-06
-1.25900106383e-06
-1.25900106383e-06
-1.25900106383e-06
-1.43887012045e-06
-1.43887012045e-06
-1.43884267629e-06
-1.43887012045e-06
8.99290394761e-07
7.19448782308e-07
5.39579725689e-07
3.59710669071e-07
1.79869056618e-07
0.0
-1.79869056618e-07
-1.79869056618e-07
-5.39579725689e-07
-7.19421338143e-07
-7.19421338143e-07
-7.19421338143e-07
-8.99290394761e-07
-1.07915945138e-06
-1.07915945138e-06
-1.07915945138e-06
-1.25900106383e-06
-1.25900106383e-06
-1.25900106383e-06
-1.43887012045e-06
-1.43887012045e-06
-1.43884267629e-06
-1.43887012045e-06
I have ported from Fortran to Python an algorithm that calculates the numerical derivative along the x direction (longitudinal) of a scalar function s on a rectilinear grid that has equal grid spacing in the x and y direction (2.5 degrees x 2.5 degrees). The input data s is oriented south to north for the latitude and west to east for the longitude. The order of the derivatives should be in the following way:
At the poles the derivatives is set to zero i.e. dsdx[j,0] = 0 and dsdx[j,-1] = 0.
Both lat and lon are one dimensional arrays, having such values:
[-90. -87.5 -85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5 -35. -32.5
-30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5 -10. -7.5 -5. -2.5
0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5 25. 27.5
30. 32.5 35. 37.5 40. 42.5 45. 47.5 50. 52.5 55. 57.5
60. 62.5 65. 67.5 70. 72.5 75. 77.5 80. 82.5 85. 87.5
90. ]
[ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5
25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5
50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5
75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5
100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5
125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5
150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5
175. 177.5 -180. -177.5 -175. -172.5 -170. -167.5 -165. -162.5
-160. -157.5 -155. -152.5 -150. -147.5 -145. -142.5 -140. -137.5
-135. -132.5 -130. -127.5 -125. -122.5 -120. -117.5 -115. -112.5
-110. -107.5 -105. -102.5 -100. -97.5 -95. -92.5 -90. -87.5
-85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5
-35. -32.5 -30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5
-10. -7.5 -5. -2.5]
8.50002
8.8
9.00002
9.20001
9.50002
9.70001
9.8
10.0
10.1
10.2
10.3
10.4
10.5
10.5
10.5
10.5
10.5
10.5
10.4
10.3
10.2
10.1
10.0
9.8
9.60001
9.40001
9.20001
9.00002
8.70001
8.50002
8.20001
7.90001
7.60001
7.3
6.90001
6.60001
8.99290394761e-07
7.19448782308e-07
5.39579725689e-07
3.59710669071e-07
1.79869056618e-07
0.0
-1.79869056618e-07
-1.79869056618e-07
-5.39579725689e-07
-7.19421338143e-07
-7.19421338143e-07
-7.19421338143e-07
-8.99290394761e-07
-1.07915945138e-06
-1.07915945138e-06
-1.07915945138e-06
-1.25900106383e-06
-1.25900106383e-06
-1.25900106383e-06
-1.43887012045e-06
-1.43887012045e-06
-1.43884267629e-06
-1.43887012045e-06
I have ported from Fortran to Python an algorithm that calculates the numerical derivative along the x direction (longitudinal) of a scalar function s on a rectilinear grid that has equal grid spacing in the x and y directions (2.5 degrees x 2.5 degrees). The input data s is oriented south to north for the latitude and west to east for the longitude. The order of the derivatives should be in the following way:
At the poles, the derivatives are set to zero, i.e., dsdx[j,0] = 0 and dsdx[j,-1] = 0.
Both lat and lon are one-dimensional arrays, having such values:
[-90. -87.5 -85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5 -35. -32.5
-30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5 -10. -7.5 -5. -2.5
0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5 25. 27.5
30. 32.5 35. 37.5 40. 42.5 45. 47.5 50. 52.5 55. 57.5
60. 62.5 65. 67.5 70. 72.5 75. 77.5 80. 82.5 85. 87.5
90. ]
[ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5
25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5
50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5
75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5
100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5
125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5
150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5
175. 177.5 -180. -177.5 -175. -172.5 -170. -167.5 -165. -162.5
-160. -157.5 -155. -152.5 -150. -147.5 -145. -142.5 -140. -137.5
-135. -132.5 -130. -127.5 -125. -122.5 -120. -117.5 -115. -112.5
-110. -107.5 -105. -102.5 -100. -97.5 -95. -92.5 -90. -87.5
-85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5
-35. -32.5 -30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5
-10. -7.5 -5. -2.5]
8.50002
8.8
9.00002
9.20001
9.50002
9.70001
9.8
10.0
10.1
10.2
10.3
10.4
10.5
10.5
10.5
10.5
10.5
10.5
10.4
10.3
10.2
10.1
10.0
9.8
9.60001
9.40001
9.20001
9.00002
8.70001
8.50002
8.20001
7.90001
7.60001
7.3
6.90001
6.60001
8.99290394761e-07
7.19448782308e-07
5.39579725689e-07
3.59710669071e-07
1.79869056618e-07
0.0
-1.79869056618e-07
-1.79869056618e-07
-5.39579725689e-07
-7.19421338143e-07
-7.19421338143e-07
-7.19421338143e-07
-8.99290394761e-07
-1.07915945138e-06
-1.07915945138e-06
-1.07915945138e-06
-1.25900106383e-06
-1.25900106383e-06
-1.25900106383e-06
-1.43887012045e-06
-1.43887012045e-06
-1.43884267629e-06
-1.43887012045e-06
Numerical differentiation on sphere with pythonPython
I have ported from Fortran to Python an algorithm that calculates the numerical derivative alongalong the x direction (longitudinal) of a scalar function ss on a rectilinear grid that has equal grid spacing in the x and y direction(2.5 degrees x 2.5 degrees). The input data ss is oriented south to north for the latitude and west to east for the longitude. The order of the derivatives should be in the following way
ds/dx = s[east] - s[west]/dlon
. At the poles the derivatives is set to zero i.e. dsdx[j,0] = 0 and dsdx[j,-1] = 0:
Both lat and lon are one dimensional arrays -
ds/dx = s[east] - s[west]/dlon
having such valuesAt the poles the derivatives is set to zero i.e. -dsdx[j,0] = 0 and dsdx[j,-1] = 0.
[-90. -87.5 -85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5 -35. -32.5
-30. -27.5 -25. -22.5 -20. -17.5Both -15. -12.5lat and -10. -7.5 -5. -2.5
0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5 25. 27.5
30. 32.5 35. 37.5 40. 42.5 45. 47.5 50. 52.5 55. 57.5
60. 62.5 65. 67.5 70. 72.5 75. 77.5 80. 82.5 85. 87.5
90. ]lon are one dimensional arrays, having such values:
and lon's values are
[-90. -87.5 -85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5 -35. -32.5
-30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5 -10. -7.5 -5. -2.5
0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5 25. 27.5
30. 32.5 35. 37.5 40. 42.5 45. 47.5 50. 52.5 55. 57.5
60. 62.5 65. 67.5 70. 72.5 75. 77.5 80. 82.5 85. 87.5
90. ]
[ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5
25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5
50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5
75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5
100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5
125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5
150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5
175. 177.5 -180. -177.5 -175. -172.5 -170. -167.5 -165. -162.5
-160. -157.5 -155. -152.5 -150. -147.5 -145. -142.5 -140. -137.5
-135. -132.5 -130. -127.5 -125. -122.5 -120. -117.5 -115. -112.5
-110. -107.5 -105. -102.5 -100. -97.5 -95. -92.5 -90. -87.5
-85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5
-35. -32.5 -30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5
-10. -7.5 -5. -2.5]lon's values are:
Where I am looking for improvement is the handling of the missing data and perhaps exception handling and obviously code clarity and any other improvements. Eventually I plan to handle the missing data via interpolation. I can share the original Fortran code for comparison as well. The constants(or globals) will eventually go to a separate class file that holds all the constants. I am including them for completion.
[ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5
25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5
50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5
75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5
100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5
125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5
150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5
175. 177.5 -180. -177.5 -175. -172.5 -170. -167.5 -165. -162.5
-160. -157.5 -155. -152.5 -150. -147.5 -145. -142.5 -140. -137.5
-135. -132.5 -130. -127.5 -125. -122.5 -120. -117.5 -115. -112.5
-110. -107.5 -105. -102.5 -100. -97.5 -95. -92.5 -90. -87.5
-85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5
-35. -32.5 -30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5
-10. -7.5 -5. -2.5]
HereWhere I am looking for improvement is some samplethe handling of the missing data and perhaps exception handling and obviously code clarity and any other improvements. Eventually I plan to handle the missing data via interpolation. I can share the original Fortran code for scomparison as well. The constants (or globals) will eventually go to a separate class file that holds all the constants. I am including them for completion.
8.50002
8.8
9.00002
9.20001
9.50002
9.70001
9.8
10.0
10.1
10.2
10.3
10.4
10.5
10.5
10.5
10.5
10.5
10.5
10.4
10.3
10.2
10.1
10.0
9.8
9.60001
9.40001
9.20001
9.00002
8.70001
8.50002
8.20001
7.90001
7.60001
7.3
6.90001
6.60001Here is some sample data for s:
and output data looks like
8.50002
8.8
9.00002
9.20001
9.50002
9.70001
9.8
10.0
10.1
10.2
10.3
10.4
10.5
10.5
10.5
10.5
10.5
10.5
10.4
10.3
10.2
10.1
10.0
9.8
9.60001
9.40001
9.20001
9.00002
8.70001
8.50002
8.20001
7.90001
7.60001
7.3
6.90001
6.60001
8.99290394761e-07 7.19448782308e-07 5.39579725689e-07 3.59710669071e-07 1.79869056618e-07 0.0 -1.79869056618e-07 -1.79869056618e-07 -5.39579725689e-07 -7.19421338143e-07 -7.19421338143e-07 -7.19421338143e-07 -8.99290394761e-07 -1.07915945138e-06 -1.07915945138e-06 -1.07915945138e-06 -1.25900106383e-06 -1.25900106383e-06 -1.25900106383e-06 -1.43887012045e-06 -1.43887012045e-06 -1.43884267629e-06 -1.43887012045e-06and output data looks like:
8.99290394761e-07
7.19448782308e-07
5.39579725689e-07
3.59710669071e-07
1.79869056618e-07
0.0
-1.79869056618e-07
-1.79869056618e-07
-5.39579725689e-07
-7.19421338143e-07
-7.19421338143e-07
-7.19421338143e-07
-8.99290394761e-07
-1.07915945138e-06
-1.07915945138e-06
-1.07915945138e-06
-1.25900106383e-06
-1.25900106383e-06
-1.25900106383e-06
-1.43887012045e-06
-1.43887012045e-06
-1.43884267629e-06
-1.43887012045e-06
Numerical differentiation on sphere with python
I have ported from Fortran to Python an algorithm that calculates the numerical derivative along the x direction (longitudinal) of a scalar function s on a rectilinear grid that has equal grid spacing in the x and y direction(2.5 degrees x 2.5 degrees). The input data s is oriented south to north for the latitude and west to east for the longitude. The order of the derivatives should be in the following way ds/dx = s[east] - s[west]/dlon . At the poles the derivatives is set to zero i.e. dsdx[j,0] = 0 and dsdx[j,-1] = 0
Both lat and lon are one dimensional arrays -
having such values -
[-90. -87.5 -85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5 -60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5 -35. -32.5 -30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5 -10. -7.5 -5. -2.5 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5 25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5 50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5 75. 77.5 80. 82.5 85. 87.5 90. ]
and lon's values are
[ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5 25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5 50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5 75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5 100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5 125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5 150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5 175. 177.5 -180. -177.5 -175. -172.5 -170. -167.5 -165. -162.5 -160. -157.5 -155. -152.5 -150. -147.5 -145. -142.5 -140. -137.5 -135. -132.5 -130. -127.5 -125. -122.5 -120. -117.5 -115. -112.5 -110. -107.5 -105. -102.5 -100. -97.5 -95. -92.5 -90. -87.5 -85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5 -60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5 -35. -32.5 -30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5 -10. -7.5 -5. -2.5]
Where I am looking for improvement is the handling of the missing data and perhaps exception handling and obviously code clarity and any other improvements. Eventually I plan to handle the missing data via interpolation. I can share the original Fortran code for comparison as well. The constants(or globals) will eventually go to a separate class file that holds all the constants. I am including them for completion.
Here is some sample data for s
8.50002 8.8 9.00002 9.20001 9.50002 9.70001 9.8 10.0 10.1 10.2 10.3 10.4 10.5 10.5 10.5 10.5 10.5 10.5 10.4 10.3 10.2 10.1 10.0 9.8 9.60001 9.40001 9.20001 9.00002 8.70001 8.50002 8.20001 7.90001 7.60001 7.3 6.90001 6.60001
and output data looks like
8.99290394761e-07 7.19448782308e-07 5.39579725689e-07 3.59710669071e-07 1.79869056618e-07 0.0 -1.79869056618e-07 -1.79869056618e-07 -5.39579725689e-07 -7.19421338143e-07 -7.19421338143e-07 -7.19421338143e-07 -8.99290394761e-07 -1.07915945138e-06 -1.07915945138e-06 -1.07915945138e-06 -1.25900106383e-06 -1.25900106383e-06 -1.25900106383e-06 -1.43887012045e-06 -1.43887012045e-06 -1.43884267629e-06 -1.43887012045e-06
Numerical differentiation on sphere with Python
I have ported from Fortran to Python an algorithm that calculates the numerical derivative along the x direction (longitudinal) of a scalar function s on a rectilinear grid that has equal grid spacing in the x and y direction(2.5 degrees x 2.5 degrees). The input data s is oriented south to north for the latitude and west to east for the longitude. The order of the derivatives should be in the following way:
ds/dx = s[east] - s[west]/dlon
At the poles the derivatives is set to zero i.e. dsdx[j,0] = 0 and dsdx[j,-1] = 0.
Both lat and lon are one dimensional arrays, having such values:
[-90. -87.5 -85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5 -35. -32.5
-30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5 -10. -7.5 -5. -2.5
0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5 25. 27.5
30. 32.5 35. 37.5 40. 42.5 45. 47.5 50. 52.5 55. 57.5
60. 62.5 65. 67.5 70. 72.5 75. 77.5 80. 82.5 85. 87.5
90. ]
lon's values are:
[ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5
25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5
50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5
75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5
100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5
125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5
150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5
175. 177.5 -180. -177.5 -175. -172.5 -170. -167.5 -165. -162.5
-160. -157.5 -155. -152.5 -150. -147.5 -145. -142.5 -140. -137.5
-135. -132.5 -130. -127.5 -125. -122.5 -120. -117.5 -115. -112.5
-110. -107.5 -105. -102.5 -100. -97.5 -95. -92.5 -90. -87.5
-85. -82.5 -80. -77.5 -75. -72.5 -70. -67.5 -65. -62.5
-60. -57.5 -55. -52.5 -50. -47.5 -45. -42.5 -40. -37.5
-35. -32.5 -30. -27.5 -25. -22.5 -20. -17.5 -15. -12.5
-10. -7.5 -5. -2.5]
Where I am looking for improvement is the handling of the missing data and perhaps exception handling and obviously code clarity and any other improvements. Eventually I plan to handle the missing data via interpolation. I can share the original Fortran code for comparison as well. The constants (or globals) will eventually go to a separate class file that holds all the constants. I am including them for completion.
Here is some sample data for s:
8.50002
8.8
9.00002
9.20001
9.50002
9.70001
9.8
10.0
10.1
10.2
10.3
10.4
10.5
10.5
10.5
10.5
10.5
10.5
10.4
10.3
10.2
10.1
10.0
9.8
9.60001
9.40001
9.20001
9.00002
8.70001
8.50002
8.20001
7.90001
7.60001
7.3
6.90001
6.60001
and output data looks like:
8.99290394761e-07
7.19448782308e-07
5.39579725689e-07
3.59710669071e-07
1.79869056618e-07
0.0
-1.79869056618e-07
-1.79869056618e-07
-5.39579725689e-07
-7.19421338143e-07
-7.19421338143e-07
-7.19421338143e-07
-8.99290394761e-07
-1.07915945138e-06
-1.07915945138e-06
-1.07915945138e-06
-1.25900106383e-06
-1.25900106383e-06
-1.25900106383e-06
-1.43887012045e-06
-1.43887012045e-06
-1.43884267629e-06
-1.43887012045e-06
import numpy as np
def ddx(s,lat,lon):
lonLen = len(lon)
latLen = len(lat)
dsdx = np.empty((latLen,lonLen))
rearth = 6371221.3
# Input data is oriented N-S. Hence invert it. It is oriented W-E.
# s is a numpy 2-D array(s(latlen,lonlen)) and I am inverting it along
# latitude axis alone -
#https://stackoverflow.com/questions/28857609/how-to-invert-the-values-of-a-two-dimensional-matrix-by-using-slicing-in-numpy
s = s[::-1,:]
# Differential increment along X direction. Equal grid spacing
# Hence is it just the difference in longitudes in radians.
# Convert into meters by multiplying the radius of earth
di = abs(np.cos(np.radians(0.0))*rearth*(np.radians(lon[1]-lon[0])))
# Missing data in grid is assigned value of -999.99
# GRID order for loop- S-N(OUTER)
# W-E(INNER)
# Loop east along a row then north
for j in range(0,latLen):
for i in range(1, lonLen-1):
if (abs(lat[j]) >= 90.0):
dsdx[j,0] = 0.0
elif (s[j, i+1] > -999.99 and s[j,i-1] > -999.99):
# ds/dx = s[east] - s[west]/2*di
dsdx[j, i] = (s[j,i+1] - s[j,i-1])/(2.*di)
elif (s[i+1,j] < -999.99 and s[j,i-1] > -999.99 and s[j,i] > -999.99):
dsdx[j,i] = (s[j,i] - s[j,i-1])/di
elif (s[j,i+1] > -999.99 and s[j,i-1] < -999.99 and s[j,i] > -999.99):
dsdx[j,i] = (s[j,i+1] - s[j,i])/di
else:
dsdx[j,i] = -999.99
# Check if the data is global and compute difference at the beginning
# and end of row J
for j in range(0,latLen):
if (abs(lat[j]) >= 90.0):
dsdx[j,0] = 0.0
dsdx[j,-1] =0.0
elif (np.allclose(2*lon[0]-lon[-1],lon[1],1e-3) or np.allclose(2*lon[0]-lon[-1],lon[1] + 360.0,1e-3)):
if (s[j, 1] > -999.99 and s[j, -1] > -999.99):
dsdx[j, 0] = (s[j, 1] - s[j,-1]) / (2.*di)
elif (s[j,1] < -999.99 and s[j,-1] > -999.99 and s[j,0] > -999.99) :
dsdx[j,0] = (s[j,1] - s[j,-1]) / di
elif (s[j, 1] > -999.99 and s[j,-1] > -999.99 and s[j,0] > -999.99):
dsdx[j,0] = (s [j,1] - s[j,0]) /di
else:
dsdx[j, 0] = -999.99
if (s[j,0] > -999.99 and s[j,-2] > -999.99):
dsdx[j,-1] = (s[j, 0] - s[j,-2]) / (2. * di)
elif (s[j,0] < -999.99 and s[j,-2] > -999.99 and s[j,-1] > -999.99):
dsdx[j,-1] = (s[j,-1] - s[j,-2]) / di
elif (s[j,0] > -999.99 and s[j,-2] < -999.99 and s[j,-1] > -999.99) :
dsdx[j,-1] = (s[j,1] - s[j,-1]) / di
else:
dsdx[j, -1] = -999.99
elif (np.allclose(lon[0],lon[-1],1e-3)):
if (s[j, 1] > -999.99 and s[j,-2] > -999.99) :
dsdx[j,0] = (s[j,1] - s[j,-2]) / (2. * di)
elif (s[j,1] < -999.99 and s[j,-2] > -999.99 and s[j,0] > -999.99) :
dsdx[j,0] = (s[j,0] - s[j,-2]) / di
elif (s[1, j] > -999.99 and s[- 2, j] < -999.99 and s[0, j] > -999.99):
dsdx[j,0] = (s[j,1] - s[j,0]) / di
else:
dsdx[j,0] = -999.99
dsdx[j,-1] = dsdx[j,0]
else:
if (s[j, 1] > -999.99 and s[j,0] > -999.99) :
dsdx[j,0] = (s[j,1] - s[j,0]) /di
else:
dsdx[j,0] = -999.99
if (s[j,-1] > -999.99 and s[j,-2] > -999.99):
dsdx[j,-1] = (s[j,-1] -s[j,-2]) /di
else:
dsdx[j,-1] = -999.99
return dsdx
import numpy as np
def ddx(s,lat,lon):
lonLen = len(lon)
latLen = len(lat)
dsdx = np.empty((latLen,lonLen))
rearth = 6371221.3
# Input data is oriented N-S. Hence invert it. It is oriented W-E.
# s is a 2-D array(s(latlen,lonlen)) and I am inverting it along
# latitude axis alone -
#https://stackoverflow.com/questions/28857609/how-to-invert-the-values-of-a-two-dimensional-matrix-by-using-slicing-in-numpy
s = s[::-1,:]
# Differential increment along X direction. Equal grid spacing
# Hence is it just the difference in longitudes in radians.
# Convert into meters by multiplying the radius of earth
di = abs(np.cos(np.radians(0.0))*rearth*(np.radians(lon[1]-lon[0])))
# GRID order for loop- S-N(OUTER)
# W-E(INNER)
# Loop east along a row then north
for j in range(0,latLen):
for i in range(1, lonLen-1):
if (abs(lat[j]) >= 90.0):
dsdx[j,0] = 0.0
elif (s[j, i+1] > -999.99 and s[j,i-1] > -999.99):
# ds/dx = s[east] - s[west]/2*di
dsdx[j, i] = (s[j,i+1] - s[j,i-1])/(2.*di)
elif (s[i+1,j] < -999.99 and s[j,i-1] > -999.99 and s[j,i] > -999.99):
dsdx[j,i] = (s[j,i] - s[j,i-1])/di
elif (s[j,i+1] > -999.99 and s[j,i-1] < -999.99 and s[j,i] > -999.99):
dsdx[j,i] = (s[j,i+1] - s[j,i])/di
else:
dsdx[j,i] = -999.99
# Check if the data is global and compute difference at the beginning
# and end of row J
for j in range(0,latLen):
if (abs(lat[j]) >= 90.0):
dsdx[j,0] = 0.0
dsdx[j,-1] =0.0
elif (np.allclose(2*lon[0]-lon[-1],lon[1],1e-3) or np.allclose(2*lon[0]-lon[-1],lon[1] + 360.0,1e-3)):
if (s[j, 1] > -999.99 and s[j, -1] > -999.99):
dsdx[j, 0] = (s[j, 1] - s[j,-1]) / (2.*di)
elif (s[j,1] < -999.99 and s[j,-1] > -999.99 and s[j,0] > -999.99) :
dsdx[j,0] = (s[j,1] - s[j,-1]) / di
elif (s[j, 1] > -999.99 and s[j,-1] > -999.99 and s[j,0] > -999.99):
dsdx[j,0] = (s [j,1] - s[j,0]) /di
else:
dsdx[j, 0] = -999.99
if (s[j,0] > -999.99 and s[j,-2] > -999.99):
dsdx[j,-1] = (s[j, 0] - s[j,-2]) / (2. * di)
elif (s[j,0] < -999.99 and s[j,-2] > -999.99 and s[j,-1] > -999.99):
dsdx[j,-1] = (s[j,-1] - s[j,-2]) / di
elif (s[j,0] > -999.99 and s[j,-2] < -999.99 and s[j,-1] > -999.99) :
dsdx[j,-1] = (s[j,1] - s[j,-1]) / di
else:
dsdx[j, -1] = -999.99
elif (np.allclose(lon[0],lon[-1],1e-3)):
if (s[j, 1] > -999.99 and s[j,-2] > -999.99) :
dsdx[j,0] = (s[j,1] - s[j,-2]) / (2. * di)
elif (s[j,1] < -999.99 and s[j,-2] > -999.99 and s[j,0] > -999.99) :
dsdx[j,0] = (s[j,0] - s[j,-2]) / di
elif (s[1, j] > -999.99 and s[- 2, j] < -999.99 and s[0, j] > -999.99):
dsdx[j,0] = (s[j,1] - s[j,0]) / di
else:
dsdx[j,0] = -999.99
dsdx[j,-1] = dsdx[j,0]
else:
if (s[j, 1] > -999.99 and s[j,0] > -999.99) :
dsdx[j,0] = (s[j,1] - s[j,0]) /di
else:
dsdx[j,0] = -999.99
if (s[j,-1] > -999.99 and s[j,-2] > -999.99):
dsdx[j,-1] = (s[j,-1] -s[j,-2]) /di
else:
dsdx[j,-1] = -999.99
return dsdx
import numpy as np
def ddx(s,lat,lon):
lonLen = len(lon)
latLen = len(lat)
dsdx = np.empty((latLen,lonLen))
rearth = 6371221.3
# Input data is oriented N-S. Hence invert it. It is oriented W-E.
# s is a numpy 2-D array(s(latlen,lonlen)) and I am inverting it along
# latitude axis alone -
#https://stackoverflow.com/questions/28857609/how-to-invert-the-values-of-a-two-dimensional-matrix-by-using-slicing-in-numpy
s = s[::-1,:]
# Differential increment along X direction. Equal grid spacing
# Hence is it just the difference in longitudes in radians.
# Convert into meters by multiplying the radius of earth
di = abs(np.cos(np.radians(0.0))*rearth*(np.radians(lon[1]-lon[0])))
# Missing data in grid is assigned value of -999.99
# GRID order for loop- S-N(OUTER)
# W-E(INNER)
# Loop east along a row then north
for j in range(0,latLen):
for i in range(1, lonLen-1):
if (abs(lat[j]) >= 90.0):
dsdx[j,0] = 0.0
elif (s[j, i+1] > -999.99 and s[j,i-1] > -999.99):
# ds/dx = s[east] - s[west]/2*di
dsdx[j, i] = (s[j,i+1] - s[j,i-1])/(2.*di)
elif (s[i+1,j] < -999.99 and s[j,i-1] > -999.99 and s[j,i] > -999.99):
dsdx[j,i] = (s[j,i] - s[j,i-1])/di
elif (s[j,i+1] > -999.99 and s[j,i-1] < -999.99 and s[j,i] > -999.99):
dsdx[j,i] = (s[j,i+1] - s[j,i])/di
else:
dsdx[j,i] = -999.99
# Check if the data is global and compute difference at the beginning
# and end of row J
for j in range(0,latLen):
if (abs(lat[j]) >= 90.0):
dsdx[j,0] = 0.0
dsdx[j,-1] =0.0
elif (np.allclose(2*lon[0]-lon[-1],lon[1],1e-3) or np.allclose(2*lon[0]-lon[-1],lon[1] + 360.0,1e-3)):
if (s[j, 1] > -999.99 and s[j, -1] > -999.99):
dsdx[j, 0] = (s[j, 1] - s[j,-1]) / (2.*di)
elif (s[j,1] < -999.99 and s[j,-1] > -999.99 and s[j,0] > -999.99) :
dsdx[j,0] = (s[j,1] - s[j,-1]) / di
elif (s[j, 1] > -999.99 and s[j,-1] > -999.99 and s[j,0] > -999.99):
dsdx[j,0] = (s [j,1] - s[j,0]) /di
else:
dsdx[j, 0] = -999.99
if (s[j,0] > -999.99 and s[j,-2] > -999.99):
dsdx[j,-1] = (s[j, 0] - s[j,-2]) / (2. * di)
elif (s[j,0] < -999.99 and s[j,-2] > -999.99 and s[j,-1] > -999.99):
dsdx[j,-1] = (s[j,-1] - s[j,-2]) / di
elif (s[j,0] > -999.99 and s[j,-2] < -999.99 and s[j,-1] > -999.99) :
dsdx[j,-1] = (s[j,1] - s[j,-1]) / di
else:
dsdx[j, -1] = -999.99
elif (np.allclose(lon[0],lon[-1],1e-3)):
if (s[j, 1] > -999.99 and s[j,-2] > -999.99) :
dsdx[j,0] = (s[j,1] - s[j,-2]) / (2. * di)
elif (s[j,1] < -999.99 and s[j,-2] > -999.99 and s[j,0] > -999.99) :
dsdx[j,0] = (s[j,0] - s[j,-2]) / di
elif (s[1, j] > -999.99 and s[- 2, j] < -999.99 and s[0, j] > -999.99):
dsdx[j,0] = (s[j,1] - s[j,0]) / di
else:
dsdx[j,0] = -999.99
dsdx[j,-1] = dsdx[j,0]
else:
if (s[j, 1] > -999.99 and s[j,0] > -999.99) :
dsdx[j,0] = (s[j,1] - s[j,0]) /di
else:
dsdx[j,0] = -999.99
if (s[j,-1] > -999.99 and s[j,-2] > -999.99):
dsdx[j,-1] = (s[j,-1] -s[j,-2]) /di
else:
dsdx[j,-1] = -999.99
return dsdx