Le portail francophone de la géomatique

crazyravioli · Tue 16 November 2010 15:30

Bonjour,
je suis en train d'ecrire un programme cense pouvoir lire plusieurs fichiers (des tuiles adjacentes) au format geotiff, et exporter dans un format "maison" la zone selectionnee sur mes tuiles concatenees. Les tuiles s'etendant parfois sur plusieurs fuseaux UTM, j'ai pris le parti de reprojeter les regions que j'extrais de chaque tuile dans la projection de la tuile "principale" (celle ou se situe le centre de la zone que j'extrais).

Je me sers donc de la methode WarpRegionToBuffer de la classe GDALWarpOperation, qui correspond a mes besoins : reprojette convenablement mon fichier d'entree et stocke le resultat dans un buffer.
Souci: WarpRegionToBuffer ne remplit pas mon buffer, celui ci reste rempli de -999.0 (valeur que jai assignee avant d'appeler WarpRegionToBuffer). La valeur de retour de WarpregionToBuffer est Ce_none, ce qui signifie que tout est cense avoir marche. de meme, CPLLasterrorMsg me renvoie une chaine vide.

Voici mon code (toutes les variables que j'utilise du genre de XSO, YSO sont initialisees dans une autre methode):

Code:

    int immai = 0;
    int jmmai = 0;    
    GDALDatasetH  hSrcDS, hDstDS;
    GDALAllRegister();
    hSrcDS = GDALOpen( this->filename.toLatin1().data(), GA_ReadOnly );
    hDstDS = GDALOpen( "out.tif", GA_Update );
    GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->hSrcDS = hSrcDS;
    psWarpOptions->hDstDS = hDstDS;
    psWarpOptions->nBandCount = 1;
    psWarpOptions->panSrcBands = 
        (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
    psWarpOptions->panSrcBands[0] = 1;
    psWarpOptions->panDstBands = 
        (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
    psWarpOptions->panDstBands[0] = 1;

    psWarpOptions->pfnProgress = GDALTermProgress;   

    // Establish reprojection transformer. 
    //    maintile->GetProjectionRef
    psWarpOptions->pTransformerArg = 
        GDALCreateGenImgProjTransformer( hSrcDS, 
        GDALGetProjectionRef(hSrcDS), 
        hDstDS,
        maintile->GetProjectionRef(), 
        FALSE, 0.0, 1 );

    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
    GDALWarpOperation oOperation;

    oOperation.Initialize( psWarpOptions );
    if (rxso < XSO)
        rxso = XSO;
    if (ryso < YSO)
        ryso = YSO;    
    if (rxne > XNE)
        rxne = XNE;
    if (ryne > YNE)
        ryne = YNE;
    immai = ABS((rxne - rxso)) / dx;
    jmmai = ABS((ryne - ryso)) / dy;
    try
    {
        this->ch2d = new float[immai * jmmai];
        if (ch2d)
        {
            for (int i = 0 ; i < immai * jmmai ; i++)
                ch2d[i] = -999;
        }
    }
    catch(...)
    {
        //could not allocate memory
        return (0);
    }


    int        nXOff=int((rxso-dx/2.-rxso)/dx);
    int     nYOff=int((YNE-(ryne+(-dy)/2.))/(dy));
    int     nXSize=((rxne-rxso+dx)/dx)+1;
    int     nYSize=((ryne-ryso+(-dy))/(dy))+1;

    //    int error =     oOperation.WarpRegionToBuffer(nXOff, nYOff, nXSize, nYSize, this->ch2d, GDT_CFloat32);
    int error =     oOperation.WarpRegionToBuffer(0, 0, nXSize, nYSize, this->ch2d,
        GDT_CFloat32, nXOff, nYOff, nXSize, nYSize);
    QString machin((char*)CPLGetLastErrorMsg());
    if (error == CE_None)
        printf("okcool\n");
    else
        printf("okpascool\n");

    GDALClose(hSrcDS);
    GDALClose(hDstDS);
    GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
    GDALDestroyWarpOptions( psWarpOptions );

    return (this->ch2d);

Mes valeurs d'offset en nombre de mailles sont correctes. Si je n'ai pas ete clair dans mes explications / code, n'hesitez pas a me le signaler...
Merci d'avance pour votre aide.

ChristopheV · Tue 16 November 2010 18:55

Bonjour,

Je ne connais pas le modèle objet GDAL, et je ne suis pas un pro du C++, mais deux pistes me viennent à l'esprit en lisant le code:
nbandcount=1 est-on sur des images noir et blanc ? (si RGB alors Nbandcount =3)
N'y a t-il pas un signe moins devant le nYsize ? (lecture des images se fait en générale de haut en bas, et (0,0) le coin sup gauche de l'image)

A+

crazyravioli · Wed 17 November 2010 10:24

Merci pour ta reponse:
J'utilise des Geotiff en noir et blanc -> nbandcount = 1
Quand j'execute mon code avec:

Code:

tuilzor.extractdomain(-64.92, 49.01, -64.5, 49.2, poDataset);

nYSize = 682, la valeur a l'air correcte (je reste dans les bornes de ma tuile).
Par contre en lisant un peu la doc de GDAL, je me suis rendu compte qu'il fallait faire un GDALWarpOperation::ChunkAndWarpImage avant de pouvoir extraire a l'aide WarpRegionToBuffer.
Voici donc mon nouveau code:

Code:

float* MyTile::extractdomain(float rxso, float ryso, float rxne, float ryne, GDALDataset *maintile)
{
    int immai = 0;
    int jmmai = 0;    
    GDALDatasetH  hSrcDS, hDstDS;
    GDALAllRegister();
    hSrcDS = GDALOpen( this->filename.toLatin1().data(), GA_ReadOnly );
    //hDstDS = GDALOpen( "out.tif", GA_Update );
    GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->hSrcDS = hSrcDS;
    psWarpOptions->hDstDS = 0;
    psWarpOptions->nBandCount = 1;
    psWarpOptions->panSrcBands = 
        (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
    psWarpOptions->panSrcBands[0] = 1;
    psWarpOptions->panDstBands = 
        (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
    psWarpOptions->panDstBands[0] = 1;

    psWarpOptions->pfnProgress = GDALTermProgress;   

    // Establish reprojection transformer. 
    //    maintile->GetProjectionRef
    psWarpOptions->pTransformerArg = 
        GDALCreateGenImgProjTransformer( hSrcDS, 
        GDALGetProjectionRef(hSrcDS), 
        0,
        maintile->GetProjectionRef(), 
        FALSE, 0.0, 1 );

    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
    GDALWarpOperation oOperation;

    oOperation.Initialize( psWarpOptions );
    if (rxso < XSO)
        rxso = XSO;
    if (ryso < YSO)
        ryso = YSO;    
    if (rxne > XNE)
        rxne = XNE;
    if (ryne > YNE)
        ryne = YNE;
    immai = ABS((rxne - rxso)) / dx;
    jmmai = ABS((ryne - ryso)) / dy;
    try
    {
        this->ch2d = new float[immai * jmmai];
        if (ch2d)
        {
            for (int i = 0 ; i < immai * jmmai ; i++)
                ch2d[i] = -999;
        }
    }
    catch(...)
    {
        //could not allocate memory
        return (0);
    }

    int        nXOff=int((rxso-dx/2.-rxso)/dx);
    int     nYOff=int((YNE-(ryne+(-dy)/2.))/(dy));
    int     nXSize=((rxne-rxso+dx)/dx)+1;
    int     nYSize=((ryne-ryso+(-dy))/(dy))+1;

    //    int error =     oOperation.WarpRegionToBuffer(nXOff, nYOff, nXSize, nYSize, this->ch2d, GDT_CFloat32);
    int xs;
    int ys;
     double adfGeoTransform[6];
    GDALSuggestedWarpOutput(hSrcDS, psWarpOptions->pfnTransformer, psWarpOptions->pTransformerArg, adfGeoTransform, &xs, &ys);
    QString machin((char*)CPLGetLastErrorMsg());

    int error2 = oOperation.ChunkAndWarpImage(0,0,xs,ys);
//    int error2 = oOperation.ChunkAndWarpImage(0,0,500,500);

    QString machin2((char*)CPLGetLastErrorMsg());

    int error =     oOperation.WarpRegionToBuffer(0, 0, nXSize, nYSize, this->ch2d,
        GDT_CFloat32, nXOff, nYOff, nXSize, nYSize);
    
    QString machin3((char*)CPLGetLastErrorMsg());
    if (error == CE_None)
        printf("okcool\n");
    else
        printf("okpascool\n");
    for (int i = 0 ; i < immai * jmmai ; i++)
    {
        if (ch2d[i] != -999)
        {
            printf("roflflflfl");
        }
    }
    GDALClose(hSrcDS);
    GDALClose(hDstDS);
    GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
    GDALDestroyWarpOptions( psWarpOptions );

    return (this->ch2d);
}

j'ai donc rajoute ceci:

Code:

GDALSuggestedWarpOutput(hSrcDS, psWarpOptions->pfnTransformer, psWarpOptions->pTransformerArg, adfGeoTransform, &xs, &ys);
    QString machin((char*)CPLGetLastErrorMsg());

    int error2 = oOperation.ChunkAndWarpImage(0,0,xs,ys);

xs et ys font bien 3601 et 3601 (taille de mon tif), par contre quand j'appelle ChunkAndWarpImage, mon appli crash : Integer division by zero

ChristopheV · Wed 17 November 2010 18:52

Bonjour,

J'ai moi aussi consulté la doc de GDAL et je me demande si pour le paramètre Width il attend la longueur en pixel ou celle en octets alignés sur un boundary de 32 bits ce qui nous ferait ici pour du noir et blanc :

(((width+7)\8) + 3) AND HFFFFFFFC

Mais ce serait plutôt une exception mémoire qui serait générée au lieu d'une division par zéro

Autre possibilité (j'ai eu un plantage similaire dans mon propre code) les header de l'image est modifié avant l'appel de la fonction et la valeur du nombre de pixel par mètre en X ou en Y est mise à Zéro.

Courage !

crazyravioli · Thu 18 November 2010 11:06

Bonjour,
d'apres ce que j'ai compris c'est bien la longueur en pixel qui est attendue.

J'ai resolu le probleme de la division par zero en compilant avec les sources de GDAL pour pouvoir y passer en mode debug, il fallait definir :

psWarpOptions->eWorkingDataType = GDT_CFloat32;
voici le code de GDAL qui crash:

Code:

 int    nWordSize = GDALGetDataTypeSize(psOptions->eWorkingDataType)/8;
[...]
if (nSrcXSize != 0 && nSrcYSize != 0 &&
        (nSrcXSize > INT_MAX / nSrcYSize ||
         nSrcXSize * nSrcYSize > INT_MAX / [b](nWordSize * psOptions->nBandCount)[/b]))
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Integer overflow : nSrcXSize=%d, nSrcYSize=%d",
                  nSrcXSize, nSrcYSize);
        return CE_Failure;
    }

si psWarpOptions->eWorkingDataType n'est pas defini, nwordsize = 0 et la c'est le drame...

toujours dans le code de WarpRregionToBuffer :

Code:

 GDALDatasetRasterIO( psOptions->hSrcDS, GF_Read, 
                                 nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize, 
                                 oWK.papabySrcImage[0], nSrcXSize, nSrcYSize,
                                 psOptions->eWorkingDataType, 
                                 psOptions->nBandCount, psOptions->panSrcBands,
                                 0, 0, 0 );

oWK.papabySrcImage[0] contient bien les altitudes de la zone que je veux extraire, mais toujours dans la projection du geotif de base.
oWK.papabyDstImage[0] est bien egal a l'adresse du pointeur que je passe en parametre, mais a aucun moment WarpRegionToBuffer ne va ecrire une seule valeur dedans, du coup le buffer dans lequel je suis cense avoir les donnees reprojetees reste rempli de la valeur par defaut (-999.0) que j'ai assignee.

Aurais-je oublie de definir quelque chose d'important dans mon instance de GDALWarpOptions ?

rouault · Sun 28 November 2010 16:11

L'utilisation du type GDT_CFloat32 me parait très douteuse. Ne serait-ce plutôt pas GDT_Float32 que tu veux utiliser ?

Petit conseil si ce n'est déjà fait : lire et repomper le code source de l'utilitaire gdalwarp ( dans apps/gdalwarp.cpp ) qui met en oeuvre les API que tu veux utiliser.

crazyravioli · Mon 29 November 2010 17:56

Salut,
j'ai pu avoir une reponse du concepteur de GDAL (frank warmerdam) via la mailing list GDAL que je vous fais partager.

----------------

Matthieu,

The core problem is that you setup a transformer that goes from pixel/line
coordinates on the source raster to georeferenced coordinates, instead of
going to pixel/line coordinates on the destination raster.

It is a clue that you use GDALSuggestedWarpOutput() to get a geotransform
and pixel size for the output buffer but you never actually make any use of
those values.

One approach is to define your own transformer instead of using
GDALCreateGenImgProjTransformer() that maps between source and
destination buffer pixel/line coordinates using the geotransforms
or perhaps even other mechanisms.

Another approach might be to create a sort of fake output file, based
on the output of GDALSuggestedWarpOutput() and then pass that in to
GDALCreateGenImgProjTransformer() but still pass NULL for the destination
dataset when you setup the psWarpOptions.

I realize using WarpRegionToBuffer() directly is rather complication - but
it was not really intended to be widely used outside of the warp engine
itself.

If you wish to followup feel free to do it by direct email to me.

Best regards,
--

Je vais donc tenter la 2eme methode, je posterai mon code en cas de reussite / probleme...

Communauté	Ressources	L'emploi	Le marché	L'association	Fils RSS
Forum Blogs Wiki Planet Sites amis	Agenda Biblio Fiches AMO-CNIG Paris PCGIS Téléchargements	Offres d'emploi Banque de CV Formations Démarche Métiers	Geo-Entreprises Geo-Communiqué Appels d'offres	L'association Qui sommes-nous ? Nous aider Bienfaiteurs Partenaires Historique Statuts	Tout le site Job Géomatique Services GeoBlogs ...	Faire un don Mentions légales

Le portail francophone de la géomatique

#1 Tue 16 November 2010 15:30

[GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage

Code:

#2 Tue 16 November 2010 18:55

Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage

#3 Wed 17 November 2010 10:24

Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage

Code:

Code:

Code:

#4 Wed 17 November 2010 18:52

Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage

#5 Thu 18 November 2010 11:06

Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage

Code:

Code:

#6 Sun 28 November 2010 16:11

Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage

#7 Mon 29 November 2010 17:56

Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage

Pied de page des forums