❶ c++,gdal,GCPS为获取地理坐标,报错求解答!
影像投影转换就是将一个地理坐标系统转换到另一个坐标系统,如果在同一个椭球基准面下的转换就是严密的转换,如果在同一个椭球体不同基准面的转换是不严密的,不同椭球体之间的转换是不严密的,这就需要用到七参数、三参数等方法。需要两个不同坐标系统下公共点坐标求得系数。例如北京54和WG4-84坐标下的同一点的经纬度或者是经过投影后的平面坐标也是不同的。那么影像投影主要分为哪些步骤呢?说白了,就三个步骤,第一,坐标转换;第二,影像的重采样,最后就是写入到新文件中。
首先来说第一步,坐标转换需要转换四个坐标,也就是四个角点。也许有人说两个点就够了,左下点和右上点。在此,我告诉你,这是错误的。投影转换后这个四个角点组成的矩形那么很有可能就不是矩形了,如果你取两个点做转换那么后面的影像投影转换后的范围就不正确了。
或者再有人问,我怎么知道影像的四个角点的坐标啊?这个很简单,通过仿射变换系数,它就是影像的像素坐标(行列号)和地理坐标之间进行关联的系数。一般是六个参数。在GDAL中可以通过以下这个函数来获得,如果影像有仿射变换系数的话。如果没有仿射变换系数但是有控制点的话也能解算出系数;如果都没有,说明这幅影像是没有地理参考的,那么狠遗憾的告诉你,这个影像就不能做投影转换。还有就是,如果你这幅影像没有投影的话也就不能做投影转换了,因为你根本不知道这幅影像的投影是啥。
CPLErr GDALDataset::GetGeoTransform
(
double *
padfTransform
)
其中padfTransform就存储了这六个参数。这是个六个double型的数的数组。
在向北向的图像中,padfTransform[1]代表像素宽度,padfTransform[5]代表像素高度。图像左上角的坐标是(padfTransform[0],padfTransform[3]),adfGeoTransform[1] X方向也就是横向的分辨率大小,padfTransform [2] 旋转系数,如果为0,就是标准的正北向图像,padfTransform [4] 旋转系数,如果为0,就是标准的正北向图像,知道了这两个参数的意义,那么我们就可以得到四个角点的地理坐标了。
在正北向的图像中,四个角点的坐标计算如下:
//计算源图像的MBR
double dbX[4];
double dbY[4];
double dbZ[4] = {0,0,0,0};
dbX[0] = adfDstGeoTransform[0]; //左上角点坐标
dbY[0] = adfDstGeoTransform[3];
//右上角坐标
dbX[1] = adfDstGeoTransform[0] + nXsize*adfDstGeoTransform[1];
dbY[1] = adfDstGeoTransform[3];
//右下角点坐标
dbX[2] = adfDstGeoTransform[0] + nXsize*adfDstGeoTransform[1] + nYsize*adfDstGeoTransform[2];
dbY[2] = adfDstGeoTransform[3] + nXsize*adfDstGeoTransform[4] + nYsize*adfDstGeoTransform[5];
//左下角坐标
dbX[3] = adfDstGeoTransform[0];
dbY[3] = adfDstGeoTransform[3] + nYsize*adfDstGeoTransform[5];
这样我们就找到了需要参与投影转换的坐标点了,下一步就是坐标转换,坐标转换的过程通过GDAL的接口实现,其底层依赖了PROJ4地图投影开源类库。对于不同的椭球体之间变换需要用到三参数布尔莎或者七参数布尔莎模型,具体过程就是首先将经纬度大地坐标转换为地心坐标系下的空间直角坐标,然后用布尔莎模型计算,最后将计算后的结果重新转换到目标地理坐标系统下的经纬度大地坐标。有了需要转换的坐标后,我们将对上述四个点的坐标进行变换,其函数如下:
[cpp] view plain
bool TranformCoordsOGR(char* pszSrcWkt,char* pszDestWkt, int nCount,double* x,double* y,double* z
,double *dfParaSrc,double* dfParaDst,int nParaCount)
{
//创建OGR的空间参考系
OGRSpatialReference oSrcSrs; //源坐标系统
OGRSpatialReference oDestSrs; //目的坐标系统
oSrcSrs.importFromWkt(&pszSrcWkt);
oDestSrs.importFromWkt(&pszDestWkt);
int nSameGeoCS = oSrcSrs.IsSameGeogCS(&oDestSrs);
//相同的椭球基准面,则进行转换
if (nSameGeoCS)
{
OGRCoordinateTransformation *poCT = NULL;
poCT = ( &oSrcSrs,&oDestSrs );
if (NULL == poCT)
{
return false;
}
int nFlag = poCT->Transform(nCount,x,y,z);
if (nFlag)
{
OGRCoordinateTransformation::DestroyCT(poCT);
return true;
}
return false;
}
else //不同的椭球体基准面,要设置七参数或者三参数
{
int nFlag = 0;
//如果是地理坐标系,直接转换为空间直角坐标
OGRErr err = 0;
double dbAsrc = 0;
double dbBsrc = 0;
double dbEsrc = 0;
dbAsrc = oSrcSrs.GetSemiMajor(&err);
dbBsrc = oSrcSrs.GetSemiMinor(&err);
dbEsrc = 1-pow((dbBsrc/dbAsrc),2.0);
if (oSrcSrs.IsProjected())
{
OGRSpatialReference* poTmpSRS = oSrcSrs.CloneGeogCS();
OGRCoordinateTransformation *poCTTmp = NULL;
poCTTmp = ( &oSrcSrs,poTmpSRS );
nFlag = poCTTmp->Transform(nCount,x,y,z);
if (!nFlag)
{
OGRCoordinateTransformation::DestroyCT(poCTTmp);
return false;
}
OGRCoordinateTransformation::DestroyCT(poCTTmp);
//pj_geodetic_to_geocentric(dbAsrc,dbEsrc,nCount,0,x,y,z);
}
//将经纬度坐标转换为空间直角坐标
CGeoEllipse geoEllipse(dbAsrc,dbBsrc);
double dbX = 0;
double dbY = 0;
double dbZ = 0;
for (int i = 0; i < nCount; i ++)
{
geoEllipse.BLH_XYZ(x[i],y[i],z[i],dbX,dbY,dbZ);
x[i] = dbX;
y[i] = dbY;
z[i] = dbZ;
}
//七参数模型
vector<double> vecX;
vector<double> vecY;
vector<double> vecZ;
for (int i = 0; i < nCount; i ++)
{
double dbX = dfParaSrc[0] + (1+dfParaSrc[6])*(x[i]+dfParaSrc[5]*y[i]-dfParaSrc[4]*z[i]);
vecX.push_back(dbX);
double dbY = dfParaSrc[1] + (1+dfParaSrc[6])*(-dfParaSrc[5]*x[i]+y[i]+dfParaSrc[3]*z[i]);
vecY.push_back(dbY);
double dbZ = dfParaSrc[2] + (1+dfParaSrc[6])*(dfParaSrc[4]*x[i]-dfParaSrc[3]*y[i]+z[i]);
vecZ.push_back(dbZ);
}
memcpy(x,&vecX[0],sizeof(double)*nCount);
memcpy(y,&vecY[0],sizeof(double)*nCount);
memcpy(z,&vecZ[0],sizeof(double)*nCount);
double dbAdst = 0;
double dbBdst = 0;
double dbEdst = 0;
dbAdst = oDestSrs.GetSemiMajor(&err);
dbBdst = oDestSrs.GetSemiMinor(&err);
//再将空间直角坐标转换为地理坐标,即经纬度坐标
CGeoEllipse geoEllipse1(dbAdst,dbBdst);
for (int i = 0; i < nCount; i ++)
{
geoEllipse1.XYZ_BLH(x[i],y[i],z[i],dbX,dbY,dbZ);
x[i] = dbX;
y[i] = dbY;
z[i] = dbZ;
}
if (oDestSrs.IsProjected())
{
const char* pszProjName = oDestSrs.GetAttrValue("PROJECTION");
OGRSpatialReference* poTmpSRS = oDestSrs.CloneGeogCS();
int nZone = oDestSrs.GetUTMZone();
char* pszTmp;
poTmpSRS->exportToWkt(&pszTmp);
OGRCoordinateTransformation *poCTTmp = NULL;
poCTTmp = ( poTmpSRS,&oDestSrs );
if (NULL == poCTTmp)
{
//MessageBox(NULL,_T("失败"),_T("提示"),MB_OK);
return false;
}
nFlag = poCTTmp->Transform(nCount,x,y,z);
if (!nFlag)
{
OGRCoordinateTransformation::DestroyCT(poCTTmp);
return false;
}
OGRCoordinateTransformation::DestroyCT(poCTTmp);
return true;
}
return true;
}
return false;
}
上述代码中有空间直角坐标和大地坐标之间的变换,这个是我自己写的,读者也可以使用PROJ中的接口进行变换。
第二步就是影像重采样了,重采样就是通过原始影像的像素值内插得到新到采样点上的像素值。这个可以直接用GDAL中重采样接口来完成。
第三步就不用详细说了,一般投影转换后需要将投影后的影像写入的新文件中,直接用GDAL的读写接口来完成。
二、GDAL影像投影转换的过程中分辨率和仿射变换系数的估算
上一节已经完成了点的投影转换,那么我们现在就要估算投影后的像素分辨率大小和仿射变换系数了。
如果投影变换前是投影坐标系统,投影转换后也是投影坐标系统,或者说另外一种情况:投影变换前是地理坐标系统,投影变换后也是地理坐标系统,并且坐标的单位都一致的,那么分辨率大小基本上没变换,可以用变换前的分辨率大小。如果变换前是地理坐标系统,投影变换后是投影坐标系统,假设地理坐标系统以度为单位,投影坐标系统以米为单位,那么投影后的像素大小可以这样估计,因为经线上一个纬度的距离大约是111km,那么变换后的分辨率可以由原始分辨率乘以111000;相反的话,如果变换前是投影坐标系统,投影变换后是地理坐标系统,假设地理坐标系统以度为单位,投影坐标系统以米为单位,同理投影后的像素大小可以这样估计,变换后的分辨率可以由原始分辨率除以111000。
仿射变换系数这样也就可以确定了,左上角的坐标就是最小x值,最大y值,在变换后的四个角点坐标中比较获得,分辨率上一段也讲了如何获得。对于正北向的图像这就够了。 然后行列数就用变换后的四个角点组成的区域的MBR的宽度除以横向分辨率得到列数,高度除以纵向分辨率得到行数。
附上出处链接:http://blog.csdn.net/zhouxuguang236/article/details/17468171
❷ 用GDAL 读取图像到Dataset,然后取出一个Band ,如何把它转换成Bitmap,便于pictureBox加载
数据 转为 byte[] ,然后用 bitmap 的 重载方法 ,有直接读取 byte[] 字节 构成 bitmap
❸ GDAL怎么用函数读取12位的图像信息
首先,采用GDAL读取图像:
GDALAllRegister();
GDALDataset*poDataset;
QStringfilename;
filename=QFileDialog::getOpenFileName(this,tr("ChooseImages"),tr("AllFles(*.*)"));
//Opentheimage
QByteArrayba=filename.toLatin1();
poDataset=(GDALDataset*)GDALOpen(ba.data(),GA_ReadOnly);
其次,成功后,可以获取图像的一些基本信息,如下:
描述信息:
const char* GDALDataset::GetDriver()->GetDescription(),通常是图像的格式;
图像大小:
图像宽度:int GDALDataset::GetRasterXSize();
图像高度:int GDALDataset::GetRasterYSize();
波段数:int GDALDataset::GetRasterCount();
投影信息:GDALDataset::GetProjectionRef(),有的图像没有投影信息;
地理坐标信息:double adfGeoTransform[6] GDALDataset::GetGeoTransform(adfGeoTransform);
波段信息:数据集中重要的信息,有波段尺寸、数据类型、颜色信息等;
获取波段的方法:poBand为指向第i个波段的指针
GDALRasterBand *poBand;
poBand = poDataset->GetRasterBand(i);
波段尺寸:
int poBand->GetXSize();
int poBand->GetYSize();
数据类型:const char* GDALGetDataTypeName(poBand->GetRasterDataType());
颜色信息:const char* (poBand->GetColorInterpretation());
❹ c语言怎么用gdal读取geotiff文件
load是导入文件,一般从mat文件中imread是图像处理工具箱的库函数,处理图像比较方便Load命令功能loadFilename将名为Filename的MAT文件中的所有变量加载到工作空间中loadFilenamexyz将名为Filename的MAT文件中的x、y、z等指定变量加载到工作空间中loadFilename-regexppat1pat2将名为Filename的MAT文件中符合表达式要求的变量加载到工作空间中loadFilenamexyz-ASCII将名为Filename的8位ASCII文件中的x、y、z等指定变量加载到工作空间中load是读取matalab本身附带的索引图(具体路径是C:\MATLAB2009\toolbox\wavelet\wavedemo);而imread是读取你自己的图片(也就是你电脑上的图)imread该函数用于读取图片文件中的数据。在matlab的命令窗口中输入docimread或者helpimread即可获得该函数的帮助信息。matlab的imread很强大,一个命令可以读取各种类型的图像。但是imread并不是一个实际功能函数。不同的图像格式有不同的编码方式,因此有不同的读取方式。实际上,为每种不同格式的图像编写各自的读取函数是适当的,实际中也是这么做的。matlab就是这样的,imread只是一个入口函数。它仅仅是做了一些文件名的处理,从你的文件名中,找到绝对路径,找到图像后缀名,然后调用合适的读取函数。比如你打开\toolbox\matlab\imagesci\private文件夹会看到很多诸如readjpg.m,readtif.m的文件。这些才是不同格式图片读取的真正函数,但是!这些函数也不是实际功能函数!你打开这些m函数就可以看到里面其实很简单。它们所做的事情和imread差不多。也是调用了一些别的函数。比如readjpg.m里的实际读取函数是rjpg8crjpg16c这些。你会发现这些文件也存在于上面所说的这个文件夹中,但是它们的后缀名不是.m,而是.mexw32(.mexw64for64bit),这些实际功能函数并不是用matlab编写的,而是用C编写的,它们是经过编译的文件,不是文本文件。matlab只是调用他们而已。也就是说实际上matlab读取图像也是调用了C语言编写的代码。而且不同格式的图像有不同的代码。imread只不过是个入口函数而已。这种结构在matlab里非常非常非常常见。管中窥豹,可见一斑,看来matlab高级语言得以应用也是建立在C语言的架构之上的
❺ VS2005平台下用GDAL做图像处理
delete poDataset; //释放资源
cout<<"RasterXSize:"<<poDataset->GetRasterXSize()<<endl;//x方向长度
cout<<"RasterYSize:"<<poDataset->GetRasterYSize()<<endl;//y方向长度
cout<<"RasterCount:"<<poDataset->GetRasterCount()<<endl;//波段数量
指针使用问题,delete语句已经把poDataset释放了,后面语句就是访问野指针了
❻ GDAL/OGR的应用
利用GDAL/OGR库,可以使基于Linux的地理空间数据管理系统提供对矢量和栅格文件数据的支持。
1 . GDAL
GDAL提供对多种栅格数据的支持,包括Arc/Info ASCII Grid(asc),GeoTiff (tiff),Erdas Imagine Images(img),ASCII DEM(dem) 等格式。1)GDAL抽象数据模型
GDAL使用抽象数据模型(abstract datamodel)来解析它所支持的数据格式,抽象数据模型包括数据集(dataset),坐标系统,仿射地理坐标转换(Affine GeoTransform), 大地控制点(GCPs), 元数据(Metadata),栅格波段(Raster Band),颜色表(ColorTable),子数据集域(Subdatasets Domain),图像结构域(Image_StructureDomain),XML域(XML:Domains)。
❼ 如何在CSharp中使用GDAL
问题解决方案,可以不考虑测试结果 将四个*_CSharp.dll在项目中“添加引用”添加进来,其余gdal16.dll和另外四个编译C#时生成的dll文件拷贝到项目的debug下。即可。 如果不把dll文件拷贝到debug下,将出现下面错误: “OSGeo.OGR.Ogr”的类型初始值设定项引发异常这样的问题。 这个问题是dll不全造成的,除了要引用的4个dll外,还有5个dll也要放到Debug目录下。 在编译C#下的gdal时,总共生成了9个dll,在编译的本机上,程序是通过环境变量path找到另外的几个dll的。 在没有编译过gdal的电脑上,反正就把这9个编译后的dll放到debug下面就一切Ok了
❽ 使用gdal库读取图像
谷歌了一下,就找到这个还好点,希望能帮到你。
http://wenku..com/view/f419727d27284b73f2425035.html
❾ 关于C++中使用头文件gdal_priv.h
你没有导入lib库,所有造成函数有定义而连接不上,将你的lib库文件加入到工程文件中,或者在原代码中加入 #pragma comment(lib,"你的lib库路径和文件名"). 然后进行编译连接。
❿ 关于UTM投影坐标与地理经纬度转换的问题
UTM是平面坐标, 转成经纬度是地理坐标