Coordinates in geostationary projection

[sfinkens]

I tried to reproduce the geographical coordinates of an MSG/SEVIRI image using pyresample.geometry.AreaDefinition.get_lonlats() but found small differences from the results obtained by the CGMS formulas for the geostationary projection. Here is an example snippet:

import numpy as np
import matplotlib.pyplot as plt
import pyresample.geometry
import pyresample.plot

# Constants
RAD2DEG = 180.0/np.pi
DEG2RAD = np.pi/180.0
EQUATOR_EARTH_RADIUS = np.float64(6378.1370)
POLE_EARTH_RADIUS = np.float64(6356.7523)
SAT_ALTITUDE = np.float64(35785.863)  # above the equator
DIST_EARTHCENTRE_SAT = EQUATOR_EARTH_RADIUS + SAT_ALTITUDE  # = 42164
COFF = 1856
LOFF = 1856
NROWS = 3712
NCOLS = 3712
CFAC = -13642337*RAD2DEG
LFAC = -13642337*RAD2DEG
ANGULAR_PIXEL_RESOLUTION = np.float64(251.53) / 3.0  # in microrad

def image2scanangle(col, row):
    """
    Convert SEVIRI image coordinates to scanning angles (also called 
    'intermediate' coordinates in the CGMS document) corresponding to the
    centre of the image pixel.

    See appendix E in
    http://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TEN_05057_SPE_MSG_LRIT_HRI&RevisionSelectionMethod=LatestReleased&Rendition=Web
    """
    x = (col - COFF) * ANGULAR_PIXEL_RESOLUTION/1E6
    y = (row - LOFF) * ANGULAR_PIXEL_RESOLUTION/1E6
    return x, y

def image2geo(col, row, lon_0=0.0):
    """
    Convert SEVIRI image coordinates to geographical lat/lon coordinates.

    See the "LRIT/HRIT Global Specification" document by CGMS, sections 4.4.3.2
    and 4.4.4:
    http://www.cgms-info.org/documents/cgms-lrit-hrit-global-specification-(v2-8-of-30-oct-2013).pdf

    @param col: SEVIRI image column(s) (Column 0 is West)
    @param row: SEVIRI image row(s) (Row 0 is North)
    @param lon_0: Nominal longitude (degrees east)
    @return: longitude(s), latitude(s)
    """
    # Compute scanning angles (or "intermediate coordinates")
    x, y = image2scanangle(col=col, row=row)

    # Transform scanning angles to geographical coordinates using the
    # geostationary projection (section 4.4.3.2)
    cosx = np.cos(x)
    cosy = np.cos(y)
    sinx = np.sin(x)
    siny = np.sin(y)

    # Some explanations to the numbers in the CGMS formulas:

    # a) The number 42164 in the formula for s_n is the distance between earth
    # centre and satellite

    # b) The number 1.006739501 in the formula for s_n and s_d is a 
    # "squared ellipsoid flatness" (EQUATOR_EARTH_RADIUS/POLE_EARTH_RADIUS)**2
    f2 = (EQUATOR_EARTH_RADIUS/POLE_EARTH_RADIUS)**2

    # c) I didn't find out where the number 1737122264 in the formula for 
    # s_d comes from...
    num = 1737122264

    sd = np.sqrt((DIST_EARTHCENTRE_SAT*cosx*cosy)**2 - (cosy**2 + f2*siny**2)*num)
    sn = (DIST_EARTHCENTRE_SAT*cosx*cosy - sd) / (cosy**2 + f2*siny**2)
    s1 = DIST_EARTHCENTRE_SAT - sn*cosx*cosy
    s2 = sn*sinx*cosy
    s3 = -sn*siny
    sxy = np.sqrt(s1**2 + s2**2)

    lon = RAD2DEG*np.arctan(s2/s1) + lon_0
    lat = RAD2DEG*np.arctan(f2*s3/sxy)

    return np.ma.masked_invalid(lon), np.ma.masked_invalid(lat)

if __name__ == '__main__':
    # Define image segment (full image here)
    col = np.arange(0, NCOLS)
    row = np.arange(0, NROWS)
    mcol, mrow = np.meshgrid(col, row)

    # Compute coordinates using CGMS formula
    lon, lat = image2geo(col=mcol, row=mrow, lon_0=0)

    # Compute coordinates using pyresample...
    #
    # ... determine projection corners. Relation between scanning angles
    # and projection coordinates is from: http://proj4.org/projections/geos.html
    x, y = image2scanangle(col=col, row=row)
    x_proj = x * SAT_ALTITUDE*1E3
    y_proj = y * SAT_ALTITUDE*1E3
    mycorners = [x_proj[0], -y_proj[-1], x_proj[-1], -y_proj[0]]  # The most eastern pixel is n=1856-1 and the most southern line is m=1856-1

    # ... setup AreaDefinition
    print mycorners
    seviri = pyresample.geometry.AreaDefinition(
        area_id='seviri',
        name='seviri',
        proj_id='seviri',
        proj_dict=dict(proj='geos',
                       lon_0=0,
                       h=SAT_ALTITUDE*1000,
                       a=EQUATOR_EARTH_RADIUS*1000,
                       b=POLE_EARTH_RADIUS*1000),
        x_size=NCOLS,
        y_size=NROWS,
        area_extent=mycorners
    )

    # ... compute coordinates and mask space pixels as they have invalid
    # values (1e30 or so)
    pyr_lon, pyr_lat = seviri.get_lonlats(dtype='float64')
    pyr_lon = np.ma.masked_where(np.logical_or(pyr_lon < -90, pyr_lon > 90),
                                 pyr_lon)
    pyr_lat = np.ma.masked_where(np.logical_or(pyr_lat < -90, pyr_lat > 90),
                                 pyr_lat)

    # Compute differences and plot results
    plt.figure()
    bmap = pyresample.plot.area_def2basemap(seviri)
    bmap.drawcoastlines()

    plt.figure()
    plt.imshow((pyr_lon - lon), cmap='bwr', vmin=-0.1, vmax=0.1)
    plt.title('lon diff')
    plt.colorbar()

    plt.figure()
    plt.imshow((pyr_lat - lat), cmap='bwr', vmin=-0.1, vmax=0.1)
    plt.title('lat diff')
    plt.colorbar()

    plt.show()

Did I understand/implement the CGMS formulas incorrectly or is that a pyresample issue? I'm using pyresample-1.5.0 .

[adybbroe]

Stephan,

I think you understood right, but the area definition in pyresample refers to the outer boundary of the area. So you have to subtract/add half a pixel on the corners.

Your mycorners are: [-5568753.6104812799, -5565753.2044411488, 5565753.2044411488, 5568753.6104812799]

But they should rather be: [-5570253.813501345, -5567253.4074612139, 5567253.4074612139, 5570253.813501345]

So your code should looks something like this:

    xscale = (x_proj[-1] - x_proj[0]) / (x_proj.shape[0] - 1)
    yscale = (y_proj[-1] - y_proj[0]) / (y_proj.shape[0] - 1)
    mycorners = [x_proj[0] - xscale / 2., -y_proj[-1] - yscale /
                 2., x_proj[-1] + xscale / 2., -y_proj[0] + yscale / 2.]

-Adam

[adybbroe]

With the above adjustment you get an area extent that differs by an offset of 5.3 meters to the one we have as default in areas.def. Here is the default definition used in pytroll:

REGION: met09globeFull {
        NAME:          Full globe MSG image 0 degrees
        PCS_ID:        geos0
        PCS_DEF:       proj=geos, lon_0=0.0, a=6378169.00, b=6356583.80, h=35785831.0
        XSIZE:         3712
        YSIZE:         3712
        AREA_EXTENT:   (-5570248.4773392612, -5567248.074173444, 5567248.074173444, 5570248.4773392612)
};

I see that a,b,h all differ slightly from yours. Not sure why!

[sfinkens]

Brilliant! Thanks for the clarification, Adam. I just copied the constants from the CGMS document for this example. In practice we use the same constants as NWCSAF-MSG does. They are slightly different from CGMS, but identical to the pytroll definition.