[gdal-dev] Re: Lost in projections

[K.-Michael Aye]

To make the sad long story short:

I was missing (=not knowing about) one method of the SpatialReference 
class: CloneGeogCS()

So if you do:

srs = osr.SpatialReference(gdaldataset.GetProjection())

and that projection is in map coordinates, to get the same thing in 
geographical coordinates one just has to do:

targetSRS = srs.CloneGeogCS()

and then to convert one into another :

transform = osr.CoordinateTransformation(srs,targetSRS)

and then for example:

lon, lat, height = transform.TransformPoint(x,y)

It is instructive to look at both srs and targetSRS.ExportToPrettyWkt() 
to see the differences.

HTH others,

Michael

On 2011-02-27 16:15:56 +0100, K.-Michael Aye said:

> 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