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
s = s[::-1,;]. Could you please correct it so that your code doesn't get closed as off-topic? \$\endgroup\$s[::-1]). More, please add your imports and some example data + what the output is. \$\endgroup\$