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