6

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:

Example image that is not an eumetsat image

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?

Cœur
  • 37,241
  • 25
  • 195
  • 267
AME
  • 2,499
  • 5
  • 29
  • 45
  • Im not familiar with pyproj but wouldn't this require "clocking" the earth based on some known locations, i.e. realizing the position of australia relative to the picture? does pyproj do this already? Or did you just mean relative lat long of a circular shape? – Paul Seeb Jun 29 '12 at 15:28
  • Sorry, I don't understand what you are asking. I want to be able to do this: for pixel in image: print latlon(pixel) – AME Jun 29 '12 at 15:41
  • 1
    How accurate do you need to be? If you can get away with a few kilometers of error, you could approximate the earth as a sphere, do some simple geometry, and forget about pyproj entirely (I don't know anything about pyproj, but I could help with the geometry). – Neil Forrester Jun 29 '12 at 15:45
  • It's for science, so as acurate as possible unfortunatelly. :/ – AME Jun 29 '12 at 15:50
  • If the satellite is geostationary, do you have a way of determining the center point of your image? – Winston Ewert Jun 29 '12 at 16:15
  • For the image, it has to be 3712/2, for the Meteosat-Satellite it is (0,0), since it is position atop of the equator and on the Greenwich prime meridian, afaik. – AME Jun 29 '12 at 16:27
  • 2
    It’s pretty easy to establish that your assumption that the satellite is above 0 longitude is wrong, because Australia is visible :-o Try again. – DisappointedByUnaccountableMod Dec 04 '18 at 08:23
  • Yes, that's an example image. Sorry for making this not more clear. But thanks for your concern. – AME Dec 05 '18 at 12:39
  • That's not a EUMETSAT image. – gerrit Jan 14 '22 at 08:49

2 Answers2

2

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.

arjenve
  • 2,313
  • 1
  • 12
  • 7
0

Firstly, this is not a EUMETSAT image. Your image shows Australia, and geostationary EUMETSAT satellites are located above (approximately) (0°N, 0°E) and (0°N, 41.5°E) for the Zero Degree Service (ZDS) and the Indian Ocean Data Coverage (IODC) service, respectively. Neither viewpoint can see Australia. Your satellite image is probably coming from the Multifunctional Transport Satellites (MTSAT) imager that was operated by the Japanese Meteorological Agency (JMA) until replaced by HIMAWARI in 2015.

To get the geolocation of pixels in images from either Meteosat SEVIRI, MTSAT, HIMAWARI AHI, or any other earth observation imager, you could also consider reading the data using the library:

from satpy import Scene
from glob import glob
filenames = glob('/path/to/your/files/*')  # for SEVIRI, can be HRIT, Native, or NetCDF
sc = Scene(filenames=filenames)
sc.load(["IR_108"], upper_right_corner="NE")
print(sc["IR_108"].attrs["area"].get_lonlat_from_array_coordinates(1000, 2000))

which gives:

(-24.259016114751823, -3.9842388279769465)

For geostationary images, this is not the fastest method, because the answer is constant, so you shouldn't need to actually read the data. However, using a well-tested and widely used library reduces the risk on error on your side. At least, it could serve as a reference. Remember that for many satellite imagers, different channels have different spatial resolutions, so the answer may be channel-dependent. The level-1 data may also be stored with the south coming before the north, such that a naïve reading and displaying shows the south at the upper part of an image and the north at the lower part, contrary to mapping conventions. Therefore, it's worth verifying that your result is accurate.

gerrit
  • 24,025
  • 17
  • 97
  • 170