[gdal-dev] Lost in projections

[K.-Michael Aye]

I just can't get my head around it... :(

I have a Mars DEM rasterdata set that is projected in polar 
stereographic with this Wkt:

PROJCS["POLAR_STEREOGRAPHIC MARS",
    GEOGCS["GCS_MARS",
        DATUM["D_MARS",
            SPHEROID["MARS",3396190,169.8944472236118]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",-90],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]

and a worldfile (output of ds.GetGeoTransform()

(-707109.6954345703,
 115.08946865797043,
 0.0,
 707109.6954345703,
 0.0,
 -115.08946865797043)

So, currently, to find the corrent sample/line coordinates to determine 
what raster window to read, I need to know the projected coordinate in 
meters and using this via gdal.ApplyGeoTransform() works fine, in both 
directions.

But I would like to find my data by providing a pair of  unprojected, 
i.e. geographical coordinates. So something like
sample, line = gdal.ApplyGeoTransform(t, lat,lon)

So I understand that I need a transformation between 
polar_stereographic and lat/long projection, but I am stumped at how I 
get a new wordfile?

I guess I can just take the above Wkt and change projection to lat/long 
and remove the parameter section, but what do I do with the affine 
transformation parameters? Where do I have to put them into the proj4 
system to achieve what I need? Somehow no examples I can find explain 
this part (or I can't see the forest because of all the trees.)

Any help would be appreciated deeply!!

Best regards,
Michael