[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