[gdal-dev] GDALSuggestedWarpOutput returns negative infinity in padfGeoTransformOut

[lancer611]

I am using a program called DS Tile that was written to take a geo referenced
image and break it into tiles to be used as input for WorldWind.  I am
modifying it slightly to fit my needs.

The problem comes when the tiling is about to begin, and a call to
GDALSuggestedWarpOutput is made:
GDALSuggestedWarpOutput(m_src, m_cSTrans->GetCSTransFunc(),
m_cSTrans->GetCSTransArg(), m_mat, &sx, &sy);

m_src is a pointer to a GDALDataset that was created with the memory driver:
GDALDataset *m_src = memDriver->Create("", xSize, ySize, 3, GDT_Byte, NULL);

The geoTransform and projection fields are set to mirror an actual
geoReferenced TIF that I have.  Of course I change the geoTransform fields 0
and 3 to reflect the offset in the original TIF.  I then use rasterio to
write the section of the TIF I want to m_src.

m_cSTrans->GetCSTransFunc() returns a pointer to a function in a class that
came with DS Tile (I have not messed with it):
int CSTrans::STransform(void *arg, int dts, int n, double *x, double *y,
double *z, int *ps) {
    
	double tx, ty;
    int i;
    CSTrans *t = reinterpret_cast<CSTrans*>(arg);
    double srcMat[6];
    double dstMat[6];
	
    if (dts) {
        for (i = 0; i < 6; ++i) srcMat[i] = t->m_dstGT.m_mat[i];
        for (i = 0; i < 6; ++i) dstMat[i] = t->m_srcIGT.m_mat[i];
    } else {
        for (i = 0; i < 6; ++i) srcMat[i] = t->m_srcGT.m_mat[i];
        for (i = 0; i < 6; ++i) dstMat[i] = t->m_dstIGT.m_mat[i];
    }
    
    for (i = 0; i < n; ++i) {
        tx = x[i]; ty = y[i];
        x[i] = srcMat[0] + tx * srcMat[1] + ty * srcMat[2];
        y[i] = srcMat[3] + tx * srcMat[4] + ty * srcMat[5];
    }

	if (t->m_ata) {
		GDALApproxTransform(t->m_ata, dts, n, x, y, z, ps);
	} else {
		GDALReprojectionTransform(t->m_rta, dts, n, x, y, z, ps);
	}

    for (i = 0; i < n; ++i) {
        tx = x[i]; ty = y[i];
        x[i] = dstMat[0] + tx * dstMat[1] + ty * dstMat[2];
        y[i] = dstMat[3] + tx * dstMat[4] + ty * dstMat[5];
    }

	return '0';  
}

m_cSTrans->GetCSTransArg() returns a pointer to an instance of the CSTrans
class (which I have not messed with):
#ifndef CSTRANS_H
#define CSTRANS_H

class CSTrans {
public:
    CSTrans();
    ~CSTrans();

public:
    void* GetCSTransArg() const { return const_cast<CSTrans*>(this); }
    GDALTransformerFunc GetCSTransFunc() const { return STransform; }

public:
    void Create(const string& srcProj, const string& dstProj, double
maxError = 0.25);
    void Destroy();

public:
    void SetSrcGT(const GeoT& gt);
    void SetDstGT(const GeoT& gt);

protected:
    static int STransform(void *arg, int dts, int n, double *x, double *y,
double *z, int *ps);

public:
    GeoT m_srcGT;
    GeoT m_srcIGT;
    GeoT m_dstGT;
    GeoT m_dstIGT;
    void *m_rta;
    void *m_ata;
};

#endif // CSTRANS_H

I think you need the create function to make sense of the STransform
function:
void CSTrans::Create(const string& srcProj, const string& dstProj, double
maxError) {
	m_rta = GDALCreateReprojectionTransformer(ProjToSRS(srcProj).c_str(),
ProjToSRS(dstProj).c_str());
	
	if (maxError >= 0.001) {
		
		m_ata = GDALCreateApproxTransformer(GDALReprojectionTransform, m_rta,
maxError);
	}
}

OK, so the problem is that the output in m_mat (double [6]) is partially
incorrect.  The second field (size in lat/lon of the tiled images) is
negative infinity (1-#IND).  This causes the tiling to not even occur
because the TIF is iterated over (somewhat) according to this field.

I know this is a lot to take it, and I will answer any questions and provide
any more code that might be needed.

Thanks in advance.

Lance
-- 
View this message in context: http://osgeo-org.1803224.n2.nabble.com/GDALSuggestedWarpOutput-returns-negative-infinity-in-padfGeoTransformOut-tp5342514p5342514.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.