1

I've downloaded EUMETSAT's RSS product (rapid-scanning-service of the northern hemisphere) - I'm interested in band 9 (IR10.8) - and I can download them as .nat (native) file or GeoTiff. Whether I download the native file and convert it into GeoTiff or use the GeoTiff as base does not change the outcome. To reduce file size I also transform into PNG/JPG which works as expected.

I'm now trying to use GDAL and PROJ to project the pixel coordinates in the image (3712 x 1392 pixel) into latitude / longitude. My reference point to validate the outcome is pixel: x, y = 1425, 666, which is the Strait of Gibraltar (quite easily detectable even from the IR imagery) which maps to a lat, lon = 35.94262669194127, -5.579847469291964

The GeoTiff contains the following projection information:

ds = gdal.Open("/folders/MSG4-SEVI-MSG15-<timestamp>.tif")
print(gdal.Info(ds))
PROJCRS["Geostationary_projection_MSG",
    BASEGEOGCRS["GCS_MSG_Ellipsoid",
        DATUM["D_MSG_DATUM",
            ELLIPSOID["MSG_SPHEROID",6356583.8,295.48806589701,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Geostationary Satellite (Sweep Y)"],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Satellite Height",35785831,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

The first odd thing I've noticed here is that in the CONVERSION section, there's a parameter:

PARAMETER["Longitude of natural origin",0, ...]

This is odd, because the RSS service should operate at 9.5 degrees east and not Greenwich.

But okay, I'm a novice at these projection definitions, so I decided to just ignore this for now. Next I tried a simple transformation using PROJ:

from pyproj import CRS, Transformer

target_crs = CRS.from_epsg(4326)  # for lat, lon
image_crs = CRS.from_wkt(ds.GetProjection())  # which should give me the `PROJCRS` above 

transformer = Transformer.from_crs(image_crs, target_crs)
lat, lon = transformer.transform(xx=1425, yy=666)
print(lat, lon)

However, this gives me (0.0060, 0.0128) (I truncated some decimal places) - this is obviously not correct.

So, I did a bit of searching and found this approach, which uses the ds.GetGeoTransform() to extract the top-left corner of the image and raster size:

(5565747.706298828, 3000.403076171875, 0.0, 1395187.4304199219, 0.0, -3000.403076171875)

The resolution (3000m / pixel) looks good, but the top-left corner doesn't seem right. So, I've loaded it into QGIS and added the OpenStreetMap layer for reference and sure enough, it seems to project the RSS image next to where Earth is supposed to be:

QGIS rendering of the GeoTiff and OpenStreetMap layer for reference

I've tried following through with this approach and end up with (inf, inf, inf). I've also tried a very similar approach based on this StackOverflow answer and got (inf, inf) as well.

Then I've noticed that the vertical centre of the RSS image seems to align pretty much with the equator and overall the image is to the right of where it should be. Given that the full-disc scan is 3712 x 3712 pixels, I thought I can just offset the pixel coordinates and "move the image" to the left and up, so it fits onto the world. So I've tried instead of passing in 1425, 666 for the X/Y pixel coordinates, I pass in 1425 - 3712, 666 - (3712 - 1392) / 2, which gives me lat/lon (27.9446, -13.6820) - at least mapping it onto Earth, but still not even close to where it should be, so the offsets are not correct like this.

I've also browsed through tons of Eumetsat documentation, but they are not written for software engineers and I could not find even a hint at why my approaches are failing. This problem seems to be very specific to EUMETSAT imagery.

1 Answer 1

0

I ended up using a method that worked reasonably well with other Geostationary satellites:

  • It relies heavily on pyproj / PROJ's built-in geostationary projection (inversed).
  • It assumes I have a full-disc image (3712 x 3712 pixels), which in my case works, since the RSS imagery shows the top 1/3 of the disc (i.e. all x/y-coordinates in the RSS image are the same as in a full-disc image).
  • I read what I can (such as satellite height and x/y resolution) from the GeoTiff and the only extra info it requires is the longitude (not included in the GeoTiff).
  • The xy_boundaries = [-5568748, 5568748] which gives ~3712 pixel at 3000m resolution.
import gdal
import numpy as np
from pyproj import CRS, Proj

# we pretend we have the full disc image (we have the top 1/3, which works for projection)
image_height, image_width = 3712, 3712 

# read resolution from dataset
ds_info = gdal.Info(gdal.Open(dataset_path), format="json")
# find out the x/y boundaries of the image (meters of the viewport from center of image)
resolution_x, resolution_y = ds_info["geoTransform"][1], ds_info["geoTransform"][5]
xy_boundaries = [round(resolution_y * image_height / 2), round(resolution_x * image_width / 2)]

# create projection
crs_info = CRS.from_wkt(ds_info["coordinateSystem"]["wkt"]).to_json_dict()
satellite_height = list(filter(lambda param: param["name"] == "Satellite Height", crs_info["conversion"]["parameters"])[0]["value"]
# geo-stationary projection with satellite height and longitude
projection = Proj(f"+proj=geos +h={satellite_height} +lon_0=9.5")

# And then create a meshgrid of x and y coordinates of the distances from the center
x, y = [], []
step_x = round((xy_boundaries[1] - xy_boundaries[0]) / image_width)
step_y = round((xy_boundaries[1] - xy_boundaries[0]) / image_height)
# generate numbers -5500000 ... +5500000 (exact number varies depending on altitude of satellite)
for i in range(xy_boundaries[0], xy_boundaries[1] + step_x, step_x):
    x.append(i)
for i in range(xy_boundaries[1], xy_boundaries[0] - step_y, -step_y):
    y.append(i)
# handle rounding issues to ensure the x and y have the correct dimension
if len(x) > image_width:
    x = x[0:image_width]
if len(y) > image_height:
    y = y[0:image_height]
# create the xx and yy matrices
xx, yy = np.meshgrid(x, y)

# create the inverse projection
lon, lat = projection(xx, yy, inverse=True)

# look-up a pixel
px_x, px_y = 1425, 666
px_lat, px_lon = lat[px_y][px_x], lon[px_y][px_x]

print("Position:", px_lat, px_lon)

Using haversine to calculate the distance to my fixed point, I end up with roughly ~1.16km distance to what I've mapped, which is not bad, given a resolution of 3km / pixel.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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