我想用GDAL写GTIFF文件,请各位高手帮忙!!

[疏影清浅]

说明：目的是想把内存中的一块数据以GTIFF的格式写出来。
*inp 为指向该块数据的指针；
*strOut 为输出的GTIFF文件名；
column 为内存中需要输出的数据块的列数；
line 为内存中需要输出的数据块的行数；
up 为内存中需要输出的数据左上角的x坐标（UTM投影）；
left 为内存中需要输出的数据左上角的y坐标（UTM投影）；
zonenumber 为内存中需要输出的数据的投影带号（UTM投影）；

图像为单波段；

现在需要把数据写成GTIFF的格式输出，UTM投影，WGS_84；
这个是我从图像中读出来的，希望程序写成的GTIFF读出来的信息和下面的相似。
string pszWKT = "PROJCS["unnamed",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",
6378137,
298.2572235630016,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",117],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1，AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32650"]]";


我刚开始学，由于马上就要交东西了，图像写不出来很着急……
下面是我写的写GTIFF图像的代码，请各位帮我看看，看看问题是什么？能帮我改通最好了，最好能给我些说明,需要的头文件
非常感谢啊！


//写GEOTIFF文件

bool WriteGEOTIFF(unsigned char *inp, char *strOut, double column,
double line, double up, double left, int zonenumber)
{
// 先注册驱动
GDALAllRegister();

// 打开数据集
GDALDataset *poDataset = (GDALDataset*)GDALOpen((LPCTSTR)inp,
GA_ReadOnly);
if( poDataset == NULL ) return false;


int nBandCount = 1;
int nXSize = line;
int nYSize = column;

// 创建GTiff格式驱动
GDALDriver *poDriver;
poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if( poDriver == NULL ) return false;

// 创建一个新数据集(这里为GTIFF影像)
char *papszParmList[] = { "INTERLEAVE = PIXEL", NULL };
GDALDataset *poDestDataset = poDriver->Create( (LPCTSTR)strOut,
nXSize, nYSize, nBandCount, GDT_Byte, papszParmList);

//设置仿射地理转换参数
double affine[6];
affine[0] = up;
affine[1] = 2.36;
affine[2] = 0;
affine[3] = left;
affine[4] = 0;
affine[5] = -2.36;

poDestDataset->SetGeoTransform(affine);


//定义一个地理坐标系
/*
OGRSpatialReference oSRS;

oSRS.SetGeogCS("My geographic coordinate system",
"WGS_1984",
"My WGS84 Spheroid",
SRS_WGS_SEMIMAJOR,
SRS_WGS_INVFLATTENING,
"Greenwich", 0.0,
"degree", SRS_UA_DEGREE_CONV);
oSRS.SetWellKnownGeogCS("WGS84");
*/

/*
string pszWKT = "PROJCS["unnamed",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",
6378137,
298.2572235630016,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",117],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1，AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32650"]]";
*/
//string pszWKT = "PROJCS["unnamed",GEOGCS["WGS 84",DATUM
["WGS_1984",SPHEROID["WGS 84",6378137,298.2572235630016,AUTHORITY
["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT
["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION
["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER
["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER
["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1，
AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32650"]]";

OGRSpatialReference oSRS;
oSRS.SetProjCS("UTM (WGS84) in northern hemisphere.");
oSRS.SetWellKnownGeogCS("WGS84");
oSRS.SetUTM(zonenumber, 1);

//poDestDataset->SetProjection(oSRS.importFromWkt(pszWKT));
//OGRSpatialReference::importFromWkt(pszWKT);

// 分配缓冲

int *lpBuf = new int[nXSize*nYSize];
//unsigned *lpBuf = new int[nXSize*nYSize];
//lpBuf = np;


// 依次读/写每个波段

//for(int i = 0; i < nBandCount; ++i)
//{
GDALRasterBand *poBand;

// 获取读影像的第i个波段

//poBand = poDataset->GetRasterBand(i+1);


// 开始位置：0, 0
// 行列数目：nXSize, nYSize

poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, lpBuf, nXSize,
nYSize, GDT_Int32, 0, 0 );

// 获取写影像的第i个波段(poDestDataset)

poBand = poDestDataset->GetRasterBand(1);

// 写第i波段的数据
// 开始位置：0, 0
// 行列数目：nXSize, nYSize

poBand->RasterIO(GF_Write, 0, 0, nXSize, nYSize, lpBuf, nXSize,
nYSize, GDT_Int32, 0, 0);
//}
// 释放缓冲

//delete[] lpBuf;

// 关闭读写文件

delete poDataset;
delete poDestDataset;

return true;
}

[li lin]

首先不要急

其次，GDALDataset *poDataset = (GDALDataset*)GDALOpen((

LPCTSTR)inp,
GA_ReadOnly);是错的，inp是内存数据，GDALOpen打开的是文件。
所以poDataset完全就是米有必要的（直接就是错的）。

所以下面poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, lpBuf, nXSize,
nYSize, GDT_Int32, 0, 0 );的read步骤也是错的。

直接在poDestDataset中用RasterIO(GF_Write就好。lpBuf换成inp就可以。

设置仿射地理转换参数没有错。

OGRSpatialReference oSRS;
oSRS.SetProjCS("UTM (WGS84) in northern hemisphere.");
oSRS.SetWellKnownGeogCS("

WGS84");
oSRS.SetUTM(zonenumber, 1);
设置也没有错。这时已经建立好了空间参考了。直接用poDestDataset->SetProjection(oSRS);就可以，不要引入wkt了。

RasterIO写数据的时候数据类型是GDT_Int32？你确定这是导师的要求？如果不是，尽量用GDT_Byte或者用GDT_UInt8

2008/12/5 疏影清浅 <air...@gmail.com>

[li lin]

头文件只要

#include "gdal_priv.h"

2008/12/5 li lin <linux23...@gmail.com>

[jun Shao]

谢谢,我先改了,调试看看

RasterIO写数据的时候数据类型是我用GDT_Byte或者用GDT_UInt8

 

2008/12/5 li lin <linux23...@gmail.com>

[jun Shao]

我按你说的改了。还有点小问题,贴出来你看看^^^,看看还有什么不对的地方?非常感谢

bool WriteGEOTIFF(unsigned char *inp, char *strOut, double column, double line, double up, double left, int zonenumber)
{

     GDALAllRegister();

     int nBandCount = 1;
     int nXSize = line;
     int nYSize = column;

    //创建GTiff格式驱动
     GDALDriver *poDriver;
     poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
     if( poDriver == NULL ) return false;

    //创建一个新数据集(这里为GTIFF影像)
    char *papszParmList[] = { "INTERLEAVE = PIXEL", NULL};

    GDALDataset *poDestDataset = poDriver->Create( (LPCTSTR)strOut,nXSize, nYSize, nBandCount, GDT_Byte, papszParmList);

   //设置仿射地理转换参数
    double affine[6];
 affine[0] = up;
 affine[1] = 2.36;
 affine[2] = 0;
 affine[3] = left;
 affine[4] = 0;
 affine[5] = -2.36;

 poDestDataset->SetGeoTransform(affine);

    //定义一个地理坐标系

   OGRSpatialReference oSRS;

   oSRS.SetProjCS("UTM (WGS84) in northern hemisphere.");
   oSRS.SetWellKnownGeogCS("WGS84");
   oSRS.SetUTM(zonenumber, 1);

   poDestDataset->SetProjection(oSRS);

   poDestDataset->RasterIO(GF_Write, 0, 0, nXSize, nYSize, inp, nXSize, nYSize, GDT_Byte, 0, 0, 0, 0, 0);

   delete poDestDataset;

    return true;
}

Compiling...
Polynomial.cpp
D:\TGRIPO\Polynomial.cpp(599) : error C2664: 'SetProjection' : cannot convert parameter 1 from 'class OGRSpatialReference' to 'const char *'
        No user-defined-conversion operator available that can perform this conversion, or the operator cannot be called
Error executing cl.exe.

Polynomial.obj - 1 error(s), 0 warning(s)

2008/12/5 jun Shao <air...@gmail.com>

[li lin]

poDestDataset->SetProjection(oSRS.ExportToWkt());

2008/12/5 jun Shao <air...@gmail.com>

[jun Shao]

谢谢,我再看看

2008/12/5 li lin <linux23...@gmail.com>

[jun Shao]

 谢谢,图象写出来了。不过WGS84,和UTM的导入还是有点问题,请高手在帮忙看看

 

另 请教,怎么给写出来的图象加上压缩啊:比如

Compression=PackBits

 

谢谢!

 

问题是:

   OGRSpatialReference oSRS;

   oSRS.SetProjCS("UTM (WGS84) in northern hemisphere.");
   oSRS.SetWellKnownGeogCS("WGS84");
   oSRS.SetUTM(zonenumber, 1);

 

poDestDataset->SetProjection(oSRS.ExportToWkt());

 

Compiling...
Polynomial.cpp
D:\TGRIPO\Polynomial.cpp(599) : error C2660: 'exportToWkt' : function does not take 0 parameters
Error executing cl.exe.

[linux23...@gmail.com]

对不起，Python写多了，
C++是这样的：

char *pszWKT = NULL;
oSRS.exportToWkt( &pszWKT );
poDestDataset->SetProjection(pszWKT);

On 12月5日, 下午2时43分, "jun Shao" <airy...@gmail.com> wrote:
> 谢谢,图象写出来了。不过WGS84,和UTM的导入还是有点问题,请高手在帮忙看看
>
> 另 请教,怎么给写出来的图象加上压缩啊:比如
> Compression=PackBits
>
> 谢谢!
>
> 问题是:
> OGRSpatialReference oSRS;
> oSRS.SetProjCS("UTM (WGS84) in northern hemisphere.");
> oSRS.SetWellKnownGeogCS("WGS84");
> oSRS.SetUTM(zonenumber, 1);
>
> poDestDataset->SetProjection(oSRS.ExportToWkt());
>
> Compiling...
> Polynomial.cpp
> D:\TGRIPO\Polynomial.cpp(599) : error C2660: 'exportToWkt' : function does
> not take 0 parameters
> Error executing cl.exe.
> Polynomial.obj - 1 error(s), 0 warning(s)
>

> 2008/12/5 jun Shao <airy...@gmail.com>
>
> > 谢谢,我再看看
>
> > 2008/12/5 li lin <linux23maill...@gmail.com>
>
> >> poDestDataset->SetProjection(oSRS.ExportToWkt());
>
> >> 2008/12/5 jun Shao <airy...@gmail.com>

> >>> 2008/12/5 jun Shao <airy...@gmail.com>
>
> >>> 谢谢,我先改了,调试看看
> >>>> RasterIO写数据的时候数据类型是我用GDT_Byte或者用GDT_UInt8
>
> >>>> 2008/12/5 li lin <linux23maill...@gmail.com>
>
> >>>>> 头文件只要
>
> >>>>> #include "gdal_priv.h"
>
> >>>>> 2008/12/5 li lin <linux23maill...@gmail.com>

>
> >>>>> 首先不要急
>
> >>>>>> 其次，GDALDataset *poDataset = (GDALDataset*)GDALOpen(( LPCTSTR)inp,
> >>>>>> GA_ReadOnly);是错的，inp是内存数据，GDALOpen打开的是文件。
> >>>>>> 所以poDataset完全就是米有必要的（直接就是错的）。
>
> >>>>>> 所以下面poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, lpBuf, nXSize,
> >>>>>> nYSize, GDT_Int32, 0, 0 );的read步骤也是错的。
>
> >>>>>> 直接在poDestDataset中用RasterIO(GF_Write就好。lpBuf换成inp就可以。
>
> >>>>>> 设置仿射地理转换参数没有错。
>
> >>>>>> OGRSpatialReference oSRS;
> >>>>>> oSRS.SetProjCS("UTM (WGS84) in northern hemisphere.");
> >>>>>> oSRS.SetWellKnownGeogCS("
> >>>>>> WGS84");
> >>>>>> oSRS.SetUTM(zonenumber, 1);
>
> >>>>>> 设置也没有错。这时已经建立好了空间参考了。直接用poDestDataset->SetProjection(oSRS);就可以，不要引入wkt了。
>
> >>>>>> RasterIO写数据的时候数据类型是GDT_Int32？你确定这是导师的要求？如果不是，尽量用GDT_Byte或者用GDT_UInt8
>

> >>>>>> 2008/12/5 疏影清浅 <airy...@gmail.com>

> >>>>>>> 我刚开始学，由于马上就要交东西了，图像写不出来很着急......

> ...
>
> 阅读更多 >>

[jun Shao]

谢谢 ！这下调试应该没有问题了，等结果出来再好好谢谢你！

对了，还有个问题是，怎么在把输出的 图像 压缩？

 

要求 采用 Compression=PackBits 方式压缩，怎么加到程序中去？

 

另，python我没用过，这个很好用吗？

2008/12/5 linux23...@gmail.com <linux23...@gmail.com>

[linux23...@gmail.com]

http://www.gdal.org/frmt_gtiff.html

看这里，有你需要的！在Creation Options里。

On 12月5日, 下午7时54分, "jun Shao" <airy...@gmail.com> wrote:
> 谢谢 ！这下调试应该没有问题了，等结果出来再好好谢谢你！
> 对了，还有个问题是，怎么在把输出的 图像 压缩？
>
> 要求 采用 Compression=PackBits 方式压缩，怎么加到程序中去？
>
> 另，python我没用过，这个很好用吗？
>

> 2008/12/5 linux23maill...@gmail.com <linux23maill...@gmail.com>

> ...
>
> 阅读更多 >>

[linux23...@gmail.com]

呵呵，至于Python，来一个，我推荐一个

On 12月5日, 下午7时54分, "jun Shao" <airy...@gmail.com> wrote:

> 谢谢 ！这下调试应该没有问题了，等结果出来再好好谢谢你！
> 对了，还有个问题是，怎么在把输出的 图像 压缩？
>
> 要求 采用 Compression=PackBits 方式压缩，怎么加到程序中去？
>
> 另，python我没用过，这个很好用吗？
>

> 2008/12/5 linux23maill...@gmail.com <linux23maill...@gmail.com>

> ...
>
> 阅读更多 >>

[jun Shao]

呵呵，好的，Python，忙过这段，好好看看，学习学习 ，呵呵，多多指教哦

2008/12/6 linux23...@gmail.com <linux23...@gmail.com>

[jun Shao]

hehe,图象写出来了,不过压缩我加的还是有问题,贴出来看看哪有问题:

 

 char *papszParmList[] = { "INTERLEAVE = PIXEL", "COMPRESS=PACKBITS", NULL};

   GDALDataset *poDestDataset = poDriver->Create( (LPCTSTR)strOut,
 nXSize, nYSize, nBandCount, GDT_Byte, papszParmList);

 

多谢拉

2008/12/6 jun Shao <air...@gmail.com>

[linux23...@gmail.com]

问题？

On 12月8日, 上午10时48分, "jun Shao" <airy...@gmail.com> wrote:
> hehe,图象写出来了,不过压缩我加的还是有问题,贴出来看看哪有问题:
>
> char *papszParmList[] = { "INTERLEAVE = PIXEL", "COMPRESS=PACKBITS", NULL};
> GDALDataset *poDestDataset = poDriver->Create( (LPCTSTR)strOut,
> nXSize, nYSize, nBandCount, GDT_Byte, papszParmList);
>
> 多谢拉
>

> 2008/12/6 jun Shao <airy...@gmail.com>

>
> > 呵呵，好的，Python，忙过这段，好好看看，学习学习 ，呵呵，多多指教哦
>

> > 2008/12/6 linux23maill...@gmail.com <linux23maill...@gmail.com>

> ...
>
> 阅读更多 >>

[jun Shao]

没有报错，只是写出来的图象没有压缩,就是设置的参数好像没有发挥作用

2008/12/8 linux23...@gmail.com <linux23...@gmail.com>

[linux23...@gmail.com]

有效果，你看不见罢了
http://www.awaresystems.be/imaging/tiff/astifftagviewer.html
这里有个tif的查看器，用这个可以看见是正确的。

On 12月8日, 下午2时07分, "jun Shao" <airy...@gmail.com> wrote:
> 没有报错，只是写出来的图象没有压缩,就是设置的参数好像没有发挥作用
>
> 2008/12/8 linux23maill...@gmail.com <linux23maill...@gmail.com>

> ...
>
> 阅读更多 >>

[jun Shao]

代码没有问题是吗?我是从文件大下看出来没有压缩的,压缩后应该小很多,不过现在和以前同样大

2008/12/8 linux23...@gmail.com <linux23...@gmail.com>

[linux23...@gmail.com]

我这里压缩了会变更大。压缩算法不见得会变小。要看数据的组织。有没有作用，看的是tif里面的标签。

On 12月9日, 上午9时11分, "jun Shao" <airy...@gmail.com> wrote:
> 代码没有问题是吗?我是从文件大下看出来没有压缩的,压缩后应该小很多,不过现在和以前同样大
>
> 2008/12/8 linux23maill...@gmail.com <linux23maill...@gmail.com>

> ...
>
> 阅读更多 >>

[jun Shao]

o ,原来这样啊，我再看看,不过若有压缩的话,用ERDAS等软件打开看属性里面应该显示的,我的显示没有压缩

 

对了,斑竹对RPC(RFM)模型  有没有 研究?用GDAL可以作网格划分和RPC 吗?

2008/12/9 linux23...@gmail.com <linux23...@gmail.com>