2

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:enter image description here

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
3
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
2
  • Yeah I understood why it was happening, I'm just curious to how other software can then compute this if this is the case when you try convert from a 2D projection to coordinates. – intern12345 Jul 25 '19 at 15:54
  • @intern12345 Can you update your question with the values you are taking from proj so your example can be reproduced? – Marcelo Villa-Piñeros Jul 25 '19 at 15:56
1

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
  ...
1
  • You are right. Grid bounds modules are greater than 0.158 for x and y in all reolutions, but the image of the Earth is an ellipse in that grid and all points that fall outside the ellipse, even if they belong to the grid, do not represent any point on Earth. Most likely, those points have null value in the grid. – Gabriel De Luca Dec 29 '19 at 5:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.