[gdal-dev] Transform from lat/long/elevation to 3d sphere

[Christopher Hunt]

Hi Peter,

Thanks for your reply - yes, that's exactly what I want to do -  
transform to cartesian coords.

Here's the math that I've already got... I was just wondering if  
GDAL's transform object somehow encapsulated it.

     latitude *= DEG_TO_RAD;
     longitude *= DEG_TO_RAD;

     double sinLatitude = sin(latitude);
     double cosLatitude = cos(latitude);

     double N = wgs84SemiMajorAxis / sqrt(1.0 - sinLatitude *  
sinLatitude * sinAngularEccentricitySq);
     double h = elevation;

     x = (N + h) * cosLatitude * cos(longitude);
     y = (N * cosAngularEccentricitySq + h) * sinLatitude;
     z = (N + h) * cosLatitude * sin(longitude);

Apparently not though given your reply. :-)

Kind regards,
Christopher

On 03/07/2009, at 5:44 PM, Peter J Halls wrote:

> Christopher,
>
>    WGS84 *is* a spherical model.  In any spherical coordinate  
> system, x & y are represented by angles to the centre of the  
> sphere.  No projection equation is going to return you x & y in  
> metres for a sphere: it makes no sense of the maths.  The 'baseline'  
> surface can be obtained from the projection specification - ie the  
> radius or ellipsoid dimensions which, with the angle, enable a  
> surface point to be calculated by simple trigonometry.  Normally,  
> one applies the local datum parameters, in order to resolve the  
> complications of the Earth's irregular surface.  The projection  
> parameters are standard and defined in the GCS table - look in the  
> GDAL/data directory for gcs.csv.
>
>    There is no projection that represents the planet in a solid, 3d  
> space using Cartesian dimensions: it simply does not make sense to  
> attempt.  Maybe it would be useful to look out the 'bible' on map  
> projections - 'Map Projections: a working manual' by the late John P  
> Snyder ... it is available for download from the USGS web site.  The  
> opening chapter of that book explains why and how projections work:  
> I think that may not only help you understand what was wrong with  
> your question but also enable you to see how to achieve your goal.   
> The later chapters give the equations for the various projections.
>
> Best wishes,
>
> Peter
>
> Christopher Hunt wrote:
>> Sorry for the hopefully not so dumb question here.
>> If I want to transform from one projection to a projection of the  
>> earth in 3D, how would I set up my target srs?
>> I want to take a WGS-84 lat/long/elevation and have an x,y,z  
>> returned in metres.
>> Here's my present code:
>>  // Read the image's projection
>>  OGRSpatialReference projSpatialReference(dataset- 
>> >GetProjectionRef());
>>  scoped_ptr<OGRSpatialReference>  
>> latLongSpatialReference(projSpatialReference.CloneGeogCS());
>>  transform =  
>> OGRCreateCoordinateTransformation(latLongSpatialReference.get(),  
>> &projSpatialReference);
>> ...and of course that transforms from lat/long/elevation to the  
>> projection srs nicely. What I'm looking for is a way to specify the  
>> target srs so that it yields the x, y, z values for WGS84 3d sphere  
>> i.e. the earth.
>> Kind regards,
>> Christopher
>> ------------------------------------------------------------------------
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- 
> --------------------------------------------------------------------------------
> Peter J Halls, GIS Advisor, University of York
> Telephone: 01904 433806     Fax: 01904 433740
> Snail mail: Computing Service, University of York, Heslington, York  
> YO10 5DD
> This message has the status of a private and personal communication
> --------------------------------------------------------------------------------

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20090703/4446598a/attachment.html