Numerical differentiation on sphere with Python
Clash Royale CLAN TAG#URR8PPP
.everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty margin-bottom:0;
up vote
1
down vote
favorite
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
python algorithm python-3.x numpy numerical-methods
edited Jan 22 at 14:22
Jamal♦
30.1k11114225
asked Jan 18 at 6:04
gansub
1186
 |Â
show 3 more comments
up vote
1
down vote
favorite
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
python algorithm python-3.x numpy numerical-methods
edited Jan 22 at 14:22
Jamal♦
30.1k11114225
asked Jan 18 at 6:04
gansub
1186
Hmm, this part doesn't look like valid Python syntax:s = s[::-1,;]. Could you please correct it so that your code doesn't get closed as off-topic?
– ÑÂүÅú
Jan 18 at 6:15
@MrGrj I made the fix. Sorry about that.
– gansub
Jan 18 at 6:17
Nope, still wrong.(PS: you might be looking for:s[::-1]). More, please add your imports and some example data + what the output is.
– ÑÂүÅú
Jan 18 at 6:19
1
The code still appears to have syntax errors, but I've retracted my close vote because after looking at the page mentioned in a comment in the code I suspect that it needs to be called with some object other than a Python list. However, the sample data makes no sense: one thing which is clear from the code is that it requires three two-dimensional arrays, and the sample data is one one-dimensional array. Maybe if you remove the list of data values and instead add some actual code which invokes the method with the appropriate types of objects it will be clearer.
– Peter Taylor
Jan 18 at 8:37
@PeterTaylor 's' is a numpy 2d array. Not a Python list. No lat and lon are one dimensional arrays. I will be happy to delete my question if there are syntax errors. I am not sure about your last sentence. Can you clarify ?
– gansub
Jan 18 at 8:46
 |Â
show 3 more comments
up vote
1
down vote
favorite
up vote
1
down vote
favorite
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
python algorithm python-3.x numpy numerical-methods
edited Jan 22 at 14:22
Jamal♦
30.1k11114225
asked Jan 18 at 6:04
gansub
1186
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
python algorithm python-3.x numpy numerical-methods
edited Jan 22 at 14:22
Jamal♦
30.1k11114225
asked Jan 18 at 6:04
gansub
1186
edited Jan 22 at 14:22
Jamal♦
30.1k11114225
edited Jan 22 at 14:22
Jamal♦
30.1k11114225
edited Jan 22 at 14:22
Jamal♦
30.1k11114225
30.1k11114225
asked Jan 18 at 6:04
gansub
1186
asked Jan 18 at 6:04
gansub
1186
asked Jan 18 at 6:04
gansub
1186
1186
Hmm, this part doesn't look like valid Python syntax:s = s[::-1,;]. Could you please correct it so that your code doesn't get closed as off-topic?
– ÑÂүÅú
Jan 18 at 6:15
@MrGrj I made the fix. Sorry about that.
– gansub
Jan 18 at 6:17
Nope, still wrong.(PS: you might be looking for:s[::-1]). More, please add your imports and some example data + what the output is.
– ÑÂүÅú
Jan 18 at 6:19
1
The code still appears to have syntax errors, but I've retracted my close vote because after looking at the page mentioned in a comment in the code I suspect that it needs to be called with some object other than a Python list. However, the sample data makes no sense: one thing which is clear from the code is that it requires three two-dimensional arrays, and the sample data is one one-dimensional array. Maybe if you remove the list of data values and instead add some actual code which invokes the method with the appropriate types of objects it will be clearer.
– Peter Taylor
Jan 18 at 8:37
@PeterTaylor 's' is a numpy 2d array. Not a Python list. No lat and lon are one dimensional arrays. I will be happy to delete my question if there are syntax errors. I am not sure about your last sentence. Can you clarify ?
– gansub
Jan 18 at 8:46
 |Â
show 3 more comments
Hmm, this part doesn't look like valid Python syntax:s = s[::-1,;]. Could you please correct it so that your code doesn't get closed as off-topic?
– ÑÂүÅú
Jan 18 at 6:15
@MrGrj I made the fix. Sorry about that.
– gansub
Jan 18 at 6:17
Nope, still wrong.(PS: you might be looking for:s[::-1]). More, please add your imports and some example data + what the output is.
– ÑÂүÅú
Jan 18 at 6:19
1
The code still appears to have syntax errors, but I've retracted my close vote because after looking at the page mentioned in a comment in the code I suspect that it needs to be called with some object other than a Python list. However, the sample data makes no sense: one thing which is clear from the code is that it requires three two-dimensional arrays, and the sample data is one one-dimensional array. Maybe if you remove the list of data values and instead add some actual code which invokes the method with the appropriate types of objects it will be clearer.
– Peter Taylor
Jan 18 at 8:37
@PeterTaylor 's' is a numpy 2d array. Not a Python list. No lat and lon are one dimensional arrays. I will be happy to delete my question if there are syntax errors. I am not sure about your last sentence. Can you clarify ?
– gansub
Jan 18 at 8:46
Hmm, this part doesn't look like valid Python syntax:
s = s[::-1,;]. Could you please correct it so that your code doesn't get closed as off-topic?– ÑÂүÅú
Jan 18 at 6:15
Hmm, this part doesn't look like valid Python syntax:
s = s[::-1,;]. Could you please correct it so that your code doesn't get closed as off-topic?– ÑÂүÅú
Jan 18 at 6:15
@MrGrj I made the fix. Sorry about that.
– gansub
Jan 18 at 6:17
@MrGrj I made the fix. Sorry about that.
– gansub
Jan 18 at 6:17
Nope, still wrong.(PS: you might be looking for:
s[::-1]). More, please add your imports and some example data + what the output is.– ÑÂүÅú
Jan 18 at 6:19
Nope, still wrong.(PS: you might be looking for:
s[::-1]). More, please add your imports and some example data + what the output is.– ÑÂүÅú
Jan 18 at 6:19
1
1
The code still appears to have syntax errors, but I've retracted my close vote because after looking at the page mentioned in a comment in the code I suspect that it needs to be called with some object other than a Python list. However, the sample data makes no sense: one thing which is clear from the code is that it requires three two-dimensional arrays, and the sample data is one one-dimensional array. Maybe if you remove the list of data values and instead add some actual code which invokes the method with the appropriate types of objects it will be clearer.
– Peter Taylor
Jan 18 at 8:37
The code still appears to have syntax errors, but I've retracted my close vote because after looking at the page mentioned in a comment in the code I suspect that it needs to be called with some object other than a Python list. However, the sample data makes no sense: one thing which is clear from the code is that it requires three two-dimensional arrays, and the sample data is one one-dimensional array. Maybe if you remove the list of data values and instead add some actual code which invokes the method with the appropriate types of objects it will be clearer.
– Peter Taylor
Jan 18 at 8:37
@PeterTaylor 's' is a numpy 2d array. Not a Python list. No lat and lon are one dimensional arrays. I will be happy to delete my question if there are syntax errors. I am not sure about your last sentence. Can you clarify ?
– gansub
Jan 18 at 8:46
@PeterTaylor 's' is a numpy 2d array. Not a Python list. No lat and lon are one dimensional arrays. I will be happy to delete my question if there are syntax errors. I am not sure about your last sentence. Can you clarify ?
– gansub
Jan 18 at 8:46
 |Â
show 3 more comments
active
oldest
votes
active
oldest
votes
active
oldest
votes
active
oldest
votes
active
oldest
votes
Â
draft saved
draft discarded
Â
draft saved
draft discarded
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fcodereview.stackexchange.com%2fquestions%2f185368%2fnumerical-differentiation-on-sphere-with-python%23new-answer', 'question_page');
);
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Hmm, this part doesn't look like valid Python syntax:
s = s[::-1,;]. Could you please correct it so that your code doesn't get closed as off-topic?– ÑÂүÅú
Jan 18 at 6:15
@MrGrj I made the fix. Sorry about that.
– gansub
Jan 18 at 6:17
Nope, still wrong.(PS: you might be looking for:
s[::-1]). More, please add your imports and some example data + what the output is.– ÑÂүÅú
Jan 18 at 6:19
1
The code still appears to have syntax errors, but I've retracted my close vote because after looking at the page mentioned in a comment in the code I suspect that it needs to be called with some object other than a Python list. However, the sample data makes no sense: one thing which is clear from the code is that it requires three two-dimensional arrays, and the sample data is one one-dimensional array. Maybe if you remove the list of data values and instead add some actual code which invokes the method with the appropriate types of objects it will be clearer.
– Peter Taylor
Jan 18 at 8:37
@PeterTaylor 's' is a numpy 2d array. Not a Python list. No lat and lon are one dimensional arrays. I will be happy to delete my question if there are syntax errors. I am not sure about your last sentence. Can you clarify ?
– gansub
Jan 18 at 8:46