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
# 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-lon[-1],lon,1e-3) or np.allclose(2*lon-lon[-1],lon + 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,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
• 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? – Grajdeanu Alex. Jan 18 '18 at 6:15
• @MrGrj I made the fix. Sorry about that. – gansub Jan 18 '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. – Grajdeanu Alex. Jan 18 '18 at 6:19
• 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 '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 '18 at 8:46