Converting GOES x,y coordinates from a 2D fixed grid into latitude and longitude

Question

I am working with some GOES data in Python and I want to grab latitude/longitude coordinates from the data. GOES' stores their projection data in radians under x,y variables that are used to project onto their fixed grid. Their PUG has a formula to convert between the two, however when I try to implement this I run in to an issue where a lot of data is lost due to a math domain error.

Here is the formula from the PUG and my code:

proj = nc.variables['goes_imager_projection']

major,minor = float(proj.semi_major_axis), float(proj.semi_minor_axis) 

H = float(proj.perspective_point_height)+major 

lambda_0 = float(proj.longitude_of_projection_origin)

a = (math.sin(x))**2 + ((math.cos(x))**2)*(((math.cos(y))**2) + 
((major**2)/minor**2)*((math.sin(y))**2))

b = -2*H*(math.cos(x))*math.cos(y)
c = H**2 - major**2

r = (-b - math.sqrt(b**2 - 4*a*c))/2*a

s_x = r*math.cos(x)*math.cos(y)
s_y = -r*math.sin(x)
s_z = r*math.cos(x)*math.sin(y)

lat = math.atan(((major**2)/minor**2)*(s_z/math.sqrt((H-s_x)**2 + s_y**2)))
lon = lambda_0 - math.atan(s_y/(H - s_x))

However, I keep getting a math domain error when it tries to compute r.

I know why it's tripping that math domain error, what I want to know is how do other programs bypass this? I know there is software that can convert from 2D grids to latitude and longitude (i.e Basemap), so what is it that my code is missing to do this?

If you wish to retry this, the values I am getting from the goes_imager_projection variable are as follows:

perspective_point_height: 35786023.0
semi_major_axis: 6378137.0
semi_minor_axis: 6356752.31414
inverse_flattening: 298.2572221
latitude_of_projection_origin: 0.0
longitude_of_projection_origin: -75.0
x = -0.106483996
y =  0.135044

Answer 1

The math domain error is probably being raised because you are passing a negative value to math.sqrt() when computing the value of r. I would check each calculation you are doing and see what b**2 - 4*a*c evaluates to because most likely it is evaluating to a negative numer.

For example, try the following code in a Python console:

import math
math.sqrt(-4)

You will see it will yield the following error:

ValueError: math domain error

Answer 2

The values of (x,y) that you are trying to evaluate are outside of the range of the Earth, hence the contents of the square root in the calculation of r is negative. These coordinates fall off in space and the corresponding projected line never intersects the Earth. With valid values of x,y, the code works well.

You could try to capture these exceptions by inserting an if statement to only execute the rest if b^2 - 4ac >= 0 :

if b**2 - 4*a*c >= 0:
  r = (-b - math.sqrt(b**2 - 4*a*c))/2*a
  ...