由经纬度计算距离和图像分辨率
如果图像只有经纬度坐标而没有投影坐标并且分辨率也是以度分秒为单位的时候,通过经纬度计算距离和图像分辨率成为很重要的问题。
我参考了该博客的文章(http://www.cnblogs.com/panyee/archive/2006/07/04/442771.html),用vc写出来是这个样子的:
1、经纬度计算两点距离:
double COilSpillOpticalDectDoc::rad(double d)
{
return d * 3.1415926 / 180.0;
}
double COilSpillOpticalDectDoc::GetDistance(double lat1, double lng1,double lat2,double lng2)
{
double radLat1 = rad(lat1);
double radLat2 = rad(lat2);
double a = radLat1 - radLat2;
double b = rad(lng1) - rad(lng2);
double s = 2 * asin(sqrt(pow(sin(a/2),2) +
cos(radLat1)*cos(radLat2)*pow(sin(b/2),2)));
s = s * 6378137;//地球半径m
s = (floor(s * 10000)) / 10000;//这个地方我不是很明白为什么要计算这一步
return s;
}
2、计算图像分辨率
GDALDataset *poDataset;
GDALAllRegister();
poDataset = (GDALDataset *) GDALOpen(env_filename, GA_ReadOnly);
nBandCount=poDataset->GetRasterCount();
int nImgSizeX=poDataset->GetRasterXSize();
int nImgSizeY=poDataset->GetRasterYSize();
double adfGeoTransform[6];
if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
{
X0=adfGeoTransform[0]; //左上角坐标
Y0=adfGeoTransform[3];
XResolution=adfGeoTransform[1]; //分辨率
YResolution=adfGeoTransform[5];
}
Raster_spf=poDataset->GetProjectionRef();//获取栅格图像的投影信息
if(Raster_spf.IsGeographic())
{
double Y_bottom=Y0+YResolution*nImgSizeY;
double Y_length=GetDistance(Y0,X0,Y_bottom,X0);
YRes=Y_length/nImgSizeY;//计算后的分辨率
double X_right=X0+XResolution*nImgSizeX;
double X_length=GetDistance(Y0,X0,Y0,X_right);
XRes=X_length/nImgSizeX;//计算后的分辨率
}
我并未进行严格的检验,不过与ArcGIS和ENVI的测量结果相比较,其计算结果还是很接近的,可以满足精度要求不是很高的需求。