[Proj] Is this correct?
Dan Crosby
dan.crosby at lincolnagritech.co.nz
Tue Jun 18 00:12:49 EST 2013
Dear List,
I'm trying to generically take the info from a GeoTIFF and get the corners in WGS84/UTM using libgeotiff and proj.
I've pieced together some code, but I don't seem to be getting the right results (negative UTM coords)
Can anyone spot any problems with this (apologies for the ugly code, I'm just hacking at the moment!):
if (m_bIsGeoTIFF)
{
x = -DBL_MAX;
y = -DBL_MAX;
double xtmp = -DBL_MAX;
double ytmp = -DBL_MAX;
switch (cnr)
{
case LowerLeft: x = 0.0; y = m_Ypx; break;
case LowerRight: x = m_Xpx; y = m_Ypx; break;
case UpperLeft: x = 0.0; y = 0.0; break;
case UpperRight: x = m_Xpx; y = 0.0; break;
}
// Try to transform the coordinate into PCS space
bRet = GTIFImageToPCS(m_pGTF, &x, &y);
if (bRet && m_GTFDef.Zone != 0)
{
xtmp = x; ytmp = y;
if (m_GTFDef.Model == ModelTypeGeographic)
{
}
else if (!GTIFProj4ToLatLong(&m_GTFDef, 1, &xtmp, &ytmp))
{
return FALSE;
}
CString spjUTM = "";
spjUTM.Format("+proj=utm +zone=%d +ellps=WGS84", UTMZone(xtmp));
projPJ pjUTM;
pjUTM = pj_init_plus(spjUTM);
if (pjUTM == NULL)
{
return FALSE;
}
char *pszProjection, **papszArgs;
projPJ psPJ;
pszProjection = GTIFGetProj4Defn(&m_GTFDef);
if (pszProjection == NULL)
{
return FALSE;
}
papszArgs = CSLTokenizeStringComplex(pszProjection, " +", TRUE, FALSE);
free(pszProjection);
psPJ = pj_init(CSLCount(papszArgs), papszArgs);
CSLDestroy(papszArgs);
if (psPJ == NULL)
{
return FALSE;
}
//x *= DEG_TO_RAD;
//y *= DEG_TO_RAD;
bRet = (pj_transform(psPJ, pjUTM, 1, 1, &x, &y, NULL) == 0);
}
}
More information about the Proj
mailing list