How to get latitude and longitude for a pixel of a eumetsat image using python?

I have a 3712x3712 pixel sized image of a geostationary eumetsat satellite. There is some black around the earth, such that the image looks like this:

For each pixel of the earth, I'd like to get its latitute and longitude. I know that there's pyproj and I was able to instantiate a projection like so:

sat = pyproj.Proj('+proj=geos +lon_0 +h=035785831.0 +x_0=0 +y_0=0')

but geting the pixel's latlon (using sat(x,y,inverse=True) where x and y are the pixel's coordinates in the image) is obviously not possible since the projection does not know the dimension (3712x3712) of my image.

What am I missing?

Solution 1:

I think you are using the correct projection library and settings.

The typical pixel resolution (km's per pixel) is reported on the eumetsat website here It reports about 3km per pixel.

You could check it by doing a lon/lat conversion to x,y on the meridian and deviding this by the pixel count (-81 degrees, 81 degrees is the maximum range, see the eumetsat site for references, http://www.eumetsat.int/):

import pyproj
sat = pyproj.Proj('+proj=geos +lon_0 +h=035785831.0 +x_0=0 +y_0=0')
x,y =  sat( 81.299, 0, radians = False, errcheck = True)
print (x * 2.0 / 3712.0 ) / 1000.0

Will give you a value of 2.927 which fits the information given by the eumetsat.

Next you can calibrate further by defining a set of well known points of your map (for example coastal features), determining their x/y position and looking up their lat/lon coordinates on-line. You could try with a range of pixel resolutions and check which one fits best or use a more elaborate routine.

The resolution may be dependent on how close you are to the equator, see here). So you may need use the above routine at several latitudes.