GDAL project TIANDITUSICHUAN tile png to wgs84 geotiff

png2tiff.py

import gdal,  osr
import numpy as np

levelsdetail = {18:(262144, 131072, 256),
                17:(131072, 65536, 256)
               }

def lat_lon2index(lat, lon, level): # top left corner
    width = levelsdetail[level][0]
    height = levelsdetail[level][1]
    # tilesize = levelsdetail[z][2]
    return int((lat+180)/(360.0/width)),int((90-lon)/(180.0/height))

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    # attention: band*rows*cols 
    cols = array.shape[2]
    rows = array.shape[1]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 3, gdal.GDT_Byte) # single band 8bit
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
 
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326) #wgs84?
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    
    outRaster.GetRasterBand(1).WriteArray(array[0,:,:][::-1])# flip it!!!
    outRaster.GetRasterBand(2).WriteArray(array[1,:,:][::-1])
    outRaster.GetRasterBand(3).WriteArray(array[2,:,:][::-1])
    outRaster.FlushCache()


def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array) # convert array to raster

if __name__ == "__main__":
    pixelWidth = 5.36441802978516E-06  # query from 
    pixelHeight = 5.36441802978516E-06 # http://www.scgis.net.cn/imap/iMapServer/defaultRest/services/newtianditudom

    lat, lon = (102.7929761, 30.17262938) # 30.12396325
    col, row= lat_lon2index(lat, lon, 18)
    rasterOrigin = (-180 + col*pixelWidth*256, 90 - (row+1)*pixelWidth*256)
    newRasterfn = '34.tif'
	
    array = gdal.Open('34.png').ReadAsArray()

    main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array)