Tile cutting Tool gdal2tiles.py rewritten as pure C + + version

gdal2tiles.py is the Python code used to generate TMS tiles in the Gdal library, supporting Google Mercator epsg:3857 and latitude and longitude epsg:4326 two kinds of tiles, output PNG format images.

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

Why should it be rewritten as pure C + +? Project Requirements Bai, a system needs to integrate a tile transduction function, but the party does not want to install complex, each time to configure the Python environment.
So began to find transduction open source resources on the Internet.
Use AE to cut graphs, call GP services directly, and use Createmapservercache, Managemapservercachetiles, and geoprocessor classes to do so. But the structure in the code is that the map service must be published first.
GeoServer in the Geowebcache middleware can also be transduction, but also need to publish the map service, and cut out the tile file naming method is very disgusting. http://www.geowebcache.org/

The core code in http://www.klokan.cz/projects/gdal2tiles/is not open source ....
In short, finally decided to rewrite gdal2tiles.py for pure C + + code.
In fact, this kind of rewriting is not complicated, gdal2tiles.py need to rewrite the code not more than 500 lines, and the Python interface called Gdal function in the C + + interface function must have, and rewrite the speed may be faster.

Here are some of the key codes that are rewritten as C + +:
According to the project needs. Only crop byte panchromatic and multi-spectral latitude and longitude projection images are supported for latitude-mesh slicing, and the initial 0-layer is two slices. Generate a JPG image.
Interface Description:

hu2tiles.exe+ + input image + + result path + + minimum number of layers + + max layers + +querysize

Where the Querysize value can take 256 or 1024, the former nearest neighbor sampling, the latter average sampling

Example:

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

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

The functions that involve coordinate transformations are as follows, and the code for Python and C + + is similar.

//////////////////////////////////////////////////////////////////////////////////////////////////
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 = (+ + lon)/Res
PY = (+ + 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::P ixelstotile (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;
}

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

Several calls to the Gdal function interface are as follows, the overall Pthon Gdal interface function is more intelligent, return to C + + a little more troublesome point ...

int Chu2tiles::hu_scale_query_to_tile (Gdaldataset *dsquery,gdaldataset *dstile)
{

int querysize = Dsquery->getrasterxsize ();
Inttilesize = 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;
}

Save the resulting image in JPG format, less than the PNG image processing an alpha band, plus not output KML files, the final C + + version of the program is faster than Python. The experimental image was reduced from 8 seconds to about 4 seconds, and more layers were not tried.

Now just rewrite the code, only to generate loose tile images, and is single-threaded processing. Subsequent considerations can be modified to multithreading.

Like this:

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

There's a: gdal2tiles-multiprocess.py.

Tile cutting Tool gdal2tiles.py rewritten as pure C + + version