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