gdal

投影系

wgs84/EPSG 4326

地心坐标系

webmercator/EPSG 3857/900913

因为这个坐标系统是 Google Map 最先使用的，或者更确切地说，是Google 最先发明的。在投影过程中，将表示地球的参考椭球体近似的作为正球体处理（正球体半径 R = 椭球体半长轴 a）。这也是为什么在 ArcGIS 中我们经常看到这个坐标系叫 WGS 1984 Web Mercator (Auxiliary Sphere)。Auxiliary Sphere 就是在告知你，这个坐标在投影过程中，将椭球体近似为正球体做投影变换，虽然基准面是WGS 1984 椭球面.

EPSG

EPSP的英文全称是European Petroleum Survey Group，中文名称为欧洲石油调查组织。这个组织成立于1986年，2005年并入IOGP(InternationalAssociation of Oil & Gas Producers)，中文名称为国际油气生产者协会.

GCJ02

GCJ-02是由中国国家测绘局（G表示Guojia国家，C表示Cehui测绘，J表示Ju局）制订的地理信息系统的坐标系统。
它是一种对经纬度数据的加密算法，即加入随机的偏差。
国内出版的各种地图系统（包括电子形式），必须至少采用GCJ-02对地理位置进行首次加密.

Bd09

百度坐标系

CGCS2000

CGCS2000与WGS84的基本定义是一致的，采用的参考椭球非常相近，椭球常数中仅扁率有细微差别，虽然因此会造成同一点在两个坐标系中的值会有微小差异，但是，在当前测量精度水平下这种微小差值是可以忽略的，因此，可以认为CGCS2000和WGS84是相容的，在坐标系的实现精度范围内两种坐标系下的坐标是一致的.

投影系转换

#coordTransform_utils.py
gcj02_to_bd09(lng, lat)
def bd09_to_gcj02(bd_lon, bd_lat)
def wgs84_to_gcj02(lng, lat)
def gcj02_to_wgs84(lng, lat)
def wgs84_to_bd09(lon, lat)
def lonlat2webmercator(lon,lat)
def webmercator2lonlat(real_x,real_y)

wgs或者gcj02坐标系的图片转为webmercator投影系

import sys
reload(sys)
sys.setdefaultencoding('utf-8') 
import gdal
import osr
import os
import math
from math import pi as PI
from coordTransform_utils import *
import cv2
# lonLat is from GPS's WGS84
# webMercator

# longitude and latitude to web Mercator
def lonlat2webmercator(lon,lat):
	x = lon
	y = math.log(math.tan((lat+90)*PI/360))/(PI/180)
	[real_x,real_y]= map(lambda x: x*20037508.34/180 , [x, y])
	return real_x,real_y

# web Mercator to longitude and latitude
def webmercator2lonlat(real_x,real_y):
	x, y = map(lambda x: x/20037508.34*180 , [real_x,real_y])
	y = 180/PI*(2*math.atan(math.exp(y*PI/180))-PI/2)
	return x, y #lon,lat

class WGS2WEBMERCATO:
	def __init__(self):
		self.driver = gdal.GetDriverByName('GTiff')
		#--- REPROJECTION ---
		self.webMercator = osr.SpatialReference()
		self.webMercator.ImportFromEPSG(3857)  # web mercator 3857
		self.wgs84 = osr.SpatialReference()
		self.wgs84.ImportFromEPSG(4326)
		
		# Work out the boundaries of the new dataset in the target projection
		self.coordTrans = osr.CoordinateTransformation(self.wgs84, self.webMercator)
		
	def wgs2webmercato(self,wgs_path,mercato_path,transform_method=1,create_method=1): 
		base_img=gdal.Open(wgs_path)
	
		xsize=base_img.RasterXSize
		ysize=base_img.RasterYSize
		band=base_img.RasterCount
	
		adfGeoTransform=base_img.GetGeoTransform()
		lon1=adfGeoTransform[0]
		lat1=adfGeoTransform[3]
		lon2=adfGeoTransform[0]+adfGeoTransform[1]*xsize
		lat2=adfGeoTransform[3]+adfGeoTransform[5]*ysize
	
		if(transform_method==1):
			
			(real_x1, real_y1, real_z1) = self.coordTrans.TransformPoint(lon1, lat1)  # corner coords in utm meters
			(real_x2, real_y2, real_z2) = self.coordTrans.TransformPoint(lon2, lat2)
			resolution_x = (real_x2 - real_x1)/(xsize) # pixel size, utm meters
			resolution_y = (real_y2 - real_y1)/(ysize)
		else:
			real_x1,real_y1=lonlat2webmercator(lon1,lat1)
			real_x2,real_y2=lonlat2webmercator(lon2,lat2)
			resolution_x=(real_x2-real_x1)/xsize
			resolution_y=(real_y2-real_y1)/ysize
	
		if(create_method==1):
			dest = self.driver.Create(mercato_path, xsize, ysize, band, gdal.GDT_Float32)
			dest.SetGeoTransform([real_x1,resolution_x,0,real_y1,0,resolution_y])
			dest.SetProjection(self.webMercator.ExportToWkt())
			gdal.ReprojectImage(base_img, dest, self.wgs84.ExportToWkt(), self.webMercator.ExportToWkt(), gdal.GRA_Bilinear)
		else:
			dest = self.driver.CreateCopy(wgs_path,base_img)
			dest.SetGeoTransform([real_x1,resolution_x,0,real_y1,0,resolution_y])
			dest.SetProjection(self.webMercator.ExportToWkt())
		
		
		for i in range(band):
			dest.GetRasterBand(i+1).SetNoDataValue(0.0)

		dest.FlushCache()  ##saves to disk!!
		
		print base_img.GetProjectionRef()
		print base_img.GetGeoTransform()
		print dest.GetProjectionRef()
		print dest.GetGeoTransform()
		
		dest = None  # Flush the dataset to the disk
		base_img = None  # only after the reprojected!
		
	def gcj022webmercator(self,base_path,mercato_path,lon1,lat1,lon2,lat2):
		base_img=gdal.Open(base_path)
	
		xsize=base_img.RasterXSize
		ysize=base_img.RasterYSize
		band=base_img.RasterCount
		
		wgsLon1,wgsLat1=gcj02_to_wgs84(lon1,lat1)
		wgsLon2,wgsLat2=gcj02_to_wgs84(lon2,lat2)
		real_x1,real_y1=lonlat2webmercator(wgsLon1,wgsLat1)
		real_x2,real_y2=lonlat2webmercator(wgsLon2,wgsLat2)
				
		resolution_x=(real_x2-real_x1)/xsize
		resolution_y=(real_y2-real_y1)/ysize
				
		dest = self.driver.CreateCopy(mercato_path,base_img)
		dest.SetProjection(self.webMercator.ExportToWkt())
		dest.SetGeoTransform([real_x1, resolution_x, 0, real_y1, 0, resolution_y])  ##sets same projection as input
		dest.FlushC