[Gdal-dev] Spaceview projection

[Enrico Zini]

Hello,

As part of the effort of geting meteosatlib[1] nearer to gdal, I have to
handle images using the spaceview projection in GDAL.  Problem is, I
didn't quote understand how.

It does not seem to be a simple affine transform.  The formula to go
from latitude,longitude to projected space that meteosatlib uses is:

Input: m.lat, m.lon.  Output: p.x, p.y:

    // sublon is the longitude of the point directly below the satellite
    // orbitRadius is the radius of the satellite orbit in units of
    //   earth radius.
    // EARTH_RPOL is the polar radius of the Earth.
    // EARTH_1E2 and EARTH_E2 are constants for which I didn't find any
    //   documentation

    // Convert to radians
    double lat = m.lat * M_PI / 180;
    double lon = (m.lon - sublon) * M_PI / 180;
    
    double c_lat;
    double r1,r2,r3,rn,rl;
    
    c_lat = atan( EARTH_1E2 * tan(lat) );
    rl    = EARTH_RPOL / ( sqrt(1.0 - EARTH_E2 * pow(cos(c_lat), 2.0)) );
    
    r1 = orbitRadius - \
          rl * cos(c_lat) * cos(lon);
    r2 = -rl * cos(c_lat) * sin(lon);
    r3 = rl * sin(c_lat);
    rn = sqrt( r1*r1 + r2*r2 + r3*r3 );
    
    p.x = atan(-r2 / r1) * 180 / M_PI;
    p.y = asin(-r3 / rn) * 180 / M_PI;

To get to the pixel coordinates in the image, a simple extra scaling
step is required, taking as input the pixel resolution and the offset of
the top left pixel:

    x = (int)rint((double)p.x * column_resolution) + image_xstart;
    y = (int)rint((double)p.y * line_resolution) + image_ystart;


Is there anything like this in GDAL?  If not, how can it be added?


Ciao,

Enrico

[1] http://meteosatlib.sourceforge.net/
-- 
GPG key: 1024D/797EBFAB 2000-12-05 Enrico Zini <enrico at debian.org>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: Digital signature
Url : http://lists.osgeo.org/pipermail/gdal-dev/attachments/20071113/4e07a1c3/attachment.bin