[Gdal-dev] RE: [vtp] converting lat/lon to terrain coords (fwd)

[dan ancona]

Hi GDALers,

I'm new on the list - the Virtual Terrain Project folks suggested that I
try asking this question here.

I'm trying to convert from some lat/lon coords coming out of our gazetteer
into the projection that the terrain I'm visualizing in VTP is using.  My
code is below, as is the text description of one of the (several) source
projections I've tried, and the target projection (below the code).

The problem is that I seem to be getting weird results out of the
tranformation. (the method returns zero)  This is a typical one...

lat 34.400830 lon -119.738330 -> -2.089828, 0.600408

Inside VTP it looks like the numbers I should be getting are in the
millions - maybe 6.00M for x and 1.97M for y.  Someone suggested that it
looks like I'm just getting the result in radians instead of feet, so I
tried setting the linear units manually (HERE in the code) even though it
looks ok in the text description.  I've tried that and a bunch of
different source projections (both WellKnown and EPSG codes) and all seem
to yield the same results.

Thanks in advance for your help.  I got this far impressively fast - C++
is not my first language and when I get something working as quickly as I
got into this, it means the documentation is good and the architecture is
clean.  Certainly seems true for GDAL.

Dan


void AClient::projectPoints(const vtProjection &targetProj){
	OGRSpatialReference sourceProj;  // will be geog
	OGRSpatialReference targetOGR;
	OGRCoordinateTransformation *trans;
	double x,y;
	DPoint3 p;
	char *ss = NULL;
	char *st = NULL;
	int res = 0;

	targetOGR = (OGRSpatialReference)targetProj;

	res = sourceProj.SetWellKnownGeogCS( "EPSG:4326" ); 
	//res = sourceProj.importFromEPSG(4269);

	// HERE trying to set linear units manually...
	res = targetOGR.SetLinearUnits(SRS_UL_US_FOOT,1.0);

    	sourceProj.exportToWkt(&ss);
	targetOGR.exportToWkt(&st);
	VTLOG("from proj: %s\nto proj:%s\n",ss,st);

	trans = OGRCreateCoordinateTransformation( &sourceProj,
						   &targetOGR);

	if(!trans)
		VTLOG("argh.");

	for (int i=0; i < m_gazList.n; ++i) {
		x = m_gazList.gazItems[i].lon;
		y = m_gazList.gazItems[i].lat;
		int result = trans->Transform(1, &x, &y);
		if (result < 1)
			VTLOG("Transformation failed.\n");
		else {
			p.x = x;	p.y = y;	p.z = 0;
			m_gazList.gazItems[i].pProj = p;
			VTLOG("lat %f lon %f -> %f, %f\n",
				m_gazList.gazItems[i].lat,m_gazList.gazItems[i].lon
				,x,y);
		}
	}
}

source projection
-----------------
GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]

target projection
-----------------
PROJCS["NAD83 / California zone 5",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS
1980",6378137,298.2572221010002,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],
UNIT["US survey foot",0.3048006096012192,
AUTHORITY["EPSG","9003"]],AUTHORITY["EPSG","26945"]]



To unsubscribe from this group, send an email to:
vtp-unsubscribe at egroups.com 

Your use of Yahoo! Groups is subject to http://docs.yahoo.com/info/terms/