瓦片切圖工具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