Numerical differentiation on sphere with Python

Question

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:

ds/dx = s[east] - s[west]/dlon

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. ]

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

• Grajdeanu Alex: Nope, still wrong.(PS: you might be looking for: s[::-1]). More, please add your imports and some example data + what the output is.
• Peter Taylor: 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.
• user74256: @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 ?
• user74256: @PeterTaylor I just ran this code with python3.5 and I got no syntax error whatsoever.
• Peter Taylor: The key word in that first clause was appears. It looks like a syntax error to anyone who knows Python but doesn't know moderately advanced features of numpy. The runtime reports it as a syntax error if it's called with standard Python lists. But if you add some code along the lines of print(ddx(None, np.create([8.50002, 8.8, 9.00002, etc]), etc)) then Python users who don't know numpy will be able to understand what's going on and review the general usage of Python. (You might even be able to do an online demo using tio.run ).

Answer 1

Documentation

The PEP 8 style guide recommends adding docstrings for functions.

The docstring should include the type of information in your question, such as:

"""
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.
"""

Naming

The PEP 8 style guide recommends snake_case for variable names. For example, lonLen would be lon_len.

rearth would be RADIUS_EARTH (upper-case since it is a constant).

Magic number

The biggest thing that stands out about the code is the repeated value -999.99. This constant should be given a meaningful name. For example:

LIMIT = -999.99

Then use it everywhere in place of the number. For example:

elif s[j, i+1] > LIMIT and s[j,i-1] > LIMIT:

Tools

You could run code development tools to automatically find some style issues with your code.

pylint shows:

C0325: Unnecessary parens after 'if' keyword (superfluous-parens)

It is common to omit parentheses for if statements. The black program can be used to automatically remove the parens.