Pages: 1
- Sujet précédent - [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage - Sujet suivant
#1 Tue 16 November 2010 15:30
- crazyravioli
- Juste Inscrit !
- Date d'inscription: 16 Nov 2010
- Messages: 4
[GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage
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.
Hors ligne
#2 Tue 16 November 2010 18:55
- ChristopheV
- Membre
- Lieu: Ajaccio
- Date d'inscription: 7 Sep 2005
- Messages: 3179
- Site web
Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage
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+
Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close
Hors ligne
#3 Wed 17 November 2010 10:24
- crazyravioli
- Juste Inscrit !
- Date d'inscription: 16 Nov 2010
- Messages: 4
Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage
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
Hors ligne
#4 Wed 17 November 2010 18:52
- ChristopheV
- Membre
- Lieu: Ajaccio
- Date d'inscription: 7 Sep 2005
- Messages: 3179
- Site web
Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage
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 !
Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close
Hors ligne
#5 Thu 18 November 2010 11:06
- crazyravioli
- Juste Inscrit !
- Date d'inscription: 16 Nov 2010
- Messages: 4
Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage
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 ?
Hors ligne
#6 Sun 28 November 2010 16:11
- rouault
- Participant assidu
- Date d'inscription: 26 Apr 2009
- Messages: 166
Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage
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.
Hors ligne
#7 Mon 29 November 2010 17:56
- crazyravioli
- Juste Inscrit !
- Date d'inscription: 16 Nov 2010
- Messages: 4
Re: [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage
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...
Hors ligne
Pages: 1
- Sujet précédent - [GDAL C++ API] utilisation de WarpRegionToBuffer pour retuilage - Sujet suivant