MyException - 我的异常网
当前位置:我的异常网» C++ » 瓦块切图工具gdal2tiles.py改写为纯c++版本

瓦块切图工具gdal2tiles.py改写为纯c++版本

www.MyException.Cn  网友分享于:2013-11-16  浏览:0次
瓦片切图工具gdal2tiles.py改写为纯c++版本

gdal2tiles.py是GDAL库中用于生成TMS瓦片的python代码,支持谷歌墨卡托EPSG:3857与经纬度EPSG:4326两种瓦片,输出png格式图像。

gdal2tiles.py More info at:

http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification
http://wiki.osgeo.org/wiki/WMS_Tiling_Client_Recommendation
http://msdn.microsoft.com/en-us/library/bb259689.aspx
http://code.google.com/apis/maps/documentation/overlays.html#Google_Maps_Coordinates

为啥要改写为纯C++的?项目需求呗,一个系统需要集成一个瓦片切图的功能,但甲方又不希望安装复杂,每次都要配置python环境。
于是开始在网上找切图的开源资源。
使用AE来切图,直接调用GP服务,利用CreateMapServerCache 、ManageMapServerCacheTiles 和Geoprocessor 类来做。但代码中的结构都是必须先发布地图服务。
GeoServer中的GeoWebCache中间件也可切图,但也是需要先发布地图服务,并且切出的瓦片文件命名方式很恶心。http://www.geowebcache.org/

http://www.klokan.cz/projects/gdal2tiles/
中核心代码不是开源的。。。。
总之最后决定改写gdal2tiles.py为纯C++代码了。
其实这种改写也不复杂,gdal2tiles.py中需要改写的代码不超过500行,并且调用的python接口gdal函数在c++接口函数里面肯定都有,并且改写后速度有可能更快。

下面是改写成C++的部分关键代码:
根据项目需要。仅支持裁切byte全色与多光谱经纬度投影图像为经纬度网格切片,初始0层为两个切片。生成jpg图像。
接口说明:

Hu2Tiles.exe+ +输入图像+ +结果路径+ +最小层数+ +最大层数+ +querysize

其中querysize数值能取256或者1024,前者最近邻采样,后者平均采样

例子:

echo %time%

Hu2Tiles.exe "D:\\GF3_MYN_QPSI_003841_E119.7_N33.2_20170503_L1A_HH_L10002340710.tiff" "D:\\huPyTiles" 2 14 256

echo %time%

pause

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

涉及到坐标转换的函数如下,可见python和C++的代码还是很相似的。

//////////////////////////////////////////////////////////////////////////////////////////////////
//def LonLatToPixels(self, lon, lat, zoom):
//"Converts lon/lat to pixel coordinates in given zoom of the EPSG:4326 pyramid"
//
// res = self.resFact / 2**zoom
// px = (180 + lon) / res
// py = (90 + lat) / res
// return px, py
void CHuGlobalGeodetic::LonLatToPixels(double lon,double lat,int zoom,double* pxy)
{
double res = resFact / pow(2,(double)zoom);
pxy[0] = (180.0 + lon) / res;
pxy[1] = (90.0 + lat) / res;
}
////////////////////////////////////////////////////////////////////////////////////////////////////
// def PixelsToTile(self, px, py):
//"Returns coordinates of the tile covering region in pixel coordinates"
//
// tx = int( math.ceil( px / float(self.tileSize) ) - 1 )
// ty = int( math.ceil( py / float(self.tileSize) ) - 1 )
// return tx, ty
void CHuGlobalGeodetic::PixelsToTile(double px,double py,int* txy)
{
txy[0] = int(ceil(px/256.0) - 1);
txy[1] = int(ceil(py/256.0) - 1);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
// def LonLatToTile(self, lon, lat, zoom):
//"Returns the tile for zoom which covers given lon/lat coordinates"
//
// px, py = self.LonLatToPixels( lon, lat, zoom)
// return self.PixelsToTile(px,py)
void CHuGlobalGeodetic::LonLatToTile(double lon,double lat,int zoom,int* txy)
{
double pxy[2] = {0.0,0.0};
double res = resFact / pow(2,(double)zoom);
pxy[0] = (180.0 + lon) / res;
pxy[1] = (90.0 + lat) / res;

txy[0] = int(ceil(pxy[0]/256.0) - 1);
txy[1] = int(ceil(pxy[1]/256.0) - 1);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
// def Resolution(self, zoom ):
//"Resolution (arc/pixel) for given zoom level (measured at Equator)"
//
// return self.resFact / 2**zoom
//#return 180 / float( 1 << (8+zoom) )
double CHuGlobalGeodetic::Resolution(int zoom)
{
return resFact / pow(2,(double)zoom);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
// def ZoomForPixelSize(self, pixelSize ):
//"Maximal scaledown zoom of the pyramid closest to the pixelSize."
//
// for i in range(MAXZOOMLEVEL):
// if pixelSize > self.Resolution(i):
// if i!=0:
// return i-1
// else:
// return 0 # We don't want to scale up
int CHuGlobalGeodetic::ZoomForPixelSize(double pixelSize)
{
for (int i=0;i<32;i++)
{
if(pixelSize > Resolution(i))
{
if (i!=0)
{
return i-1;
}
else
{
return 0;
}
}
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////
// def TileBounds(self, tx, ty, zoom):
//"Returns bounds of the given tile"
// res = self.resFact / 2**zoom
// return (
// tx*self.tileSize*res - 180,
// ty*self.tileSize*res - 90,
// (tx+1)*self.tileSize*res - 180,
// (ty+1)*self.tileSize*res - 90
// )
void CHuGlobalGeodetic::TileBounds(int tx, int ty,int zoom, double* bound4)
{
double res = resFact / pow(2,(double)zoom);
bound4[0] = tx * 256.0 * res - 180.0;
bound4[1] = ty * 256.0 * res - 90.0;
bound4[2] = (tx+1) * 256.0 * res - 180.0;
bound4[3] = (ty+1) * 256.0 * res - 90.0;
}

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

几个调用gdal函数接口的例子如下,总体上pthon的gdal接口函数更智能,换回C++的稍微麻烦点。。。

int CHu2Tiles::hu_scale_query_to_tile(GDALDataset *dsquery,GDALDataset *dstile)
{

int querysize = dsquery->GetRasterXSize();
int tilesize = dstile->GetRasterXSize();
int tilebands = dstile->GetRasterCount();

if (resampling == "average")
{
if (tilebands == 1)
{
GDALRasterBandH *pRasterBand = new GDALRasterBandH();
pRasterBand[0] = dstile->GetRasterBand(1);
GDALRegenerateOverviews(dsquery->GetRasterBand(1),1,pRasterBand,"AVERAGE",NULL,NULL);
//dstile->GetRasterBand(2)->SetNoDataValue(0);
}
if (tilebands == 3)
{
GDALRasterBandH *pRasterBand = new GDALRasterBandH();
pRasterBand[0] = dstile->GetRasterBand(1);
pRasterBand[1] = dstile->GetRasterBand(2);
pRasterBand[2] = dstile->GetRasterBand(3);
GDALRegenerateOverviews(dsquery->GetRasterBand(1),3,pRasterBand,"AVERAGE",NULL,NULL);
//dstile->GetRasterBand(4)->SetNoDataValue(0);
}
}
else
{
double trans1[6] ={0.0,tilesize/(float)querysize,0.0,0.0,0.0,tilesize/(float)querysize};
double trans2[6] ={0.0,1.0,0.0,0.0,0.0,1.0};
dsquery->SetGeoTransform(trans1);
dstile->SetGeoTransform(trans2);
GDALReprojectImage(dsquery,NULL,dstile,NULL,GRA_Bilinear,0,0,NULL,NULL,NULL);
}

return 0;
}

保存结果图像为jpg格式,就比png图像少处理了一个alpha波段,加上不输出KML文件,最终C++版本程序要比python的快些。实验图像从8秒缩减到4秒左右,更多分层的还没试。

目前只是改写代码,只能生成松散的瓦片图像,并且是单线程处理。后续可考虑修改为多线程。

比如这个:

https://github.com/commenthol/gdal2tiles-leaflet

里面有个:gdal2tiles-multiprocess.py







 
 
1楼Macy
最近在研究栅格图像,也看到了py版本的。但是不满足需求,请问能发一份c++版本的源码到我邮箱吗?感激不尽 plummoon@gmail.com

文章评论

程序员都该阅读的书
程序员都该阅读的书
如何成为一名黑客
如何成为一名黑客
那些争议最大的编程观点
那些争议最大的编程观点
初级 vs 高级开发者 哪个性价比更高?
初级 vs 高级开发者 哪个性价比更高?
Google伦敦新总部 犹如星级庄园
Google伦敦新总部 犹如星级庄园
 程序员的样子
程序员的样子
Web开发人员为什么越来越懒了?
Web开发人员为什么越来越懒了?
程序员周末都喜欢做什么?
程序员周末都喜欢做什么?
Java程序员必看电影
Java程序员必看电影
亲爱的项目经理,我恨你
亲爱的项目经理,我恨你
“肮脏的”IT工作排行榜
“肮脏的”IT工作排行榜
程序员最害怕的5件事 你中招了吗?
程序员最害怕的5件事 你中招了吗?
看13位CEO、创始人和高管如何提高工作效率
看13位CEO、创始人和高管如何提高工作效率
一个程序员的时间管理
一个程序员的时间管理
中美印日四国程序员比较
中美印日四国程序员比较
程序员的鄙视链
程序员的鄙视链
为啥Android手机总会越用越慢?
为啥Android手机总会越用越慢?
10个帮程序员减压放松的网站
10个帮程序员减压放松的网站
Web开发者需具备的8个好习惯
Web开发者需具备的8个好习惯
Java 与 .NET 的平台发展之争
Java 与 .NET 的平台发展之争
为什么程序员都是夜猫子
为什么程序员都是夜猫子
程序员应该关注的一些事儿
程序员应该关注的一些事儿
要嫁就嫁程序猿—钱多话少死的早
要嫁就嫁程序猿—钱多话少死的早
团队中“技术大拿”并非越多越好
团队中“技术大拿”并非越多越好
60个开发者不容错过的免费资源库
60个开发者不容错过的免费资源库
老程序员的下场
老程序员的下场
“懒”出效率是程序员的美德
“懒”出效率是程序员的美德
编程语言是女人
编程语言是女人
程序猿的崛起——Growth Hacker
程序猿的崛起——Growth Hacker
程序员的一天:一寸光阴一寸金
程序员的一天:一寸光阴一寸金
老美怎么看待阿里赴美上市
老美怎么看待阿里赴美上市
如何区分一个程序员是“老手“还是“新手“?
如何区分一个程序员是“老手“还是“新手“?
我跳槽是因为他们的显示器更大
我跳槽是因为他们的显示器更大
什么才是优秀的用户界面设计
什么才是优秀的用户界面设计
代码女神横空出世
代码女神横空出世
十大编程算法助程序员走上高手之路
十大编程算法助程序员走上高手之路
那些性感的让人尖叫的程序员
那些性感的让人尖叫的程序员
2013年美国开发者薪资调查报告
2013年美国开发者薪资调查报告
聊聊HTTPS和SSL/TLS协议
聊聊HTTPS和SSL/TLS协议
旅行,写作,编程
旅行,写作,编程
程序员眼里IE浏览器是什么样的
程序员眼里IE浏览器是什么样的
10个调试和排错的小建议
10个调试和排错的小建议
总结2014中国互联网十大段子
总结2014中国互联网十大段子
程序员和编码员之间的区别
程序员和编码员之间的区别
2013年中国软件开发者薪资调查报告
2013年中国软件开发者薪资调查报告
当下全球最炙手可热的八位少年创业者
当下全球最炙手可热的八位少年创业者
漫画:程序员的工作
漫画:程序员的工作
做程序猿的老婆应该注意的一些事情
做程序猿的老婆应该注意的一些事情
不懂技术不要对懂技术的人说这很容易实现
不懂技术不要对懂技术的人说这很容易实现
软件开发程序错误异常ExceptionCopyright © 2009-2015 MyException 版权所有