I am reprojecting rasters in python using GDAL. I need to project several tiffs from geographic WGS 84 coordinates to WGS 1984 Web Mercator (Auxiliary Sphere), in order to use them later in Openlayers together with OpenStreetMap and maybe Google maps. I am using Python 2.7.5 and GDAL 1.10.1 from [here](http://www.lfd.uci.edu/~gohlke/pythonlibs/), and transforming coordinates using advices from [here](http://jgomezdans.github.io/gdal_notes/reprojection.html) (my code is below). In short, I imported *osgeo.osr* and used *ImportFromEPSG(code)* and *CoordinateTransformation(from,to)*.

I tried firstly **EPSG(32629)** which is UTM zone 29 and got this projected raster (more or less fine), so the code seems to be correct:
![utm](https://i.stack.imgur.com/wVeiD.png)
Then I used **EPSG(3857)** because I've read [this](https://gis.stackexchange.com/questions/50347/which-is-the-correct-projection-in-arcgis-for-epsg3857-900913) and [this](https://gis.stackexchange.com/questions/40538/what-is-the-difference-between-epsg900913-and-epsg3857) questions and found that it is the correct recent **valid code**. But the raster is created with no spatial reference at all. It is far away up in WGS 84 data frame (but will be ok if I switch the data frame to Web Mercator).
![3857](https://i.stack.imgur.com/EQ1yM.png)

With **EPSG(900913)** the output is georeferenced, but shifted about 3 raster cells to the north:
![900913](https://i.stack.imgur.com/LUM60.png)

When I reproject the raster using **ArcGIS** (export in WGS\_1984\_Web\_Mercator\_Auxiliary\_Sphere) the result is nearly fine:
![arcgis](https://i.stack.imgur.com/zbiCE.png)

And when I use the old code **102113 (41001,54004)** the result is perfect:
![54004](https://i.stack.imgur.com/5Dk7N.png)

The summary of my tests using [all codes](http://wiki.openstreetmap.org/wiki/EPSG:3857):

```
3857: far away up (missing georeference)
3785: far away up (like 3857)
3587: far away right
900913: slightly jumped up
102100: python error
102113: perfect
41001: perfect
54004: perfect
ArcGIS (web merc. aux.): good

```

**So my questions are:** 

* Why does the correct EPSG code give me wrong results?
* And why the old codes work fine, aren't they deprecated?
* Maybe my GDAL version is not good or I have errors in my python code?

The code:

```
    yres = round(lons[1]-lons[0], 4)  # pixel size, degrees
    xres = round(lats[1]-lats[0], 4)
    ysize = len(lats)-1  # number of pixels
    xsize = len(lons)-1
    ulx = round(lons[0], 4)
    uly = round(lats[-1], 4)  # last
    driver = gdal.GetDriverByName(fileformat)
    ds = driver.Create(filename, xsize, ysize, 2, gdal.GDT_Float32)  # 2 bands
    #--- Geographic ---
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # Geographic lat/lon WGS 84
    ds.SetProjection(srs.ExportToWkt())
    gt = [ulx, xres, 0, uly, 0, -yres]  # the affine transformation coeffs (ulx, pixel, angle(skew), uly, angle, -pixel)
    ds.SetGeoTransform(gt)  # coords of top left corner of top left pixel (w-file - center of the pixel!)
    outband = ds.GetRasterBand(1)
    outband.WriteArray(data)
    outband2 = ds.GetRasterBand(2)
    outband2.WriteArray(data3)
    #--- REPROJECTION ---
    utm29 = osr.SpatialReference()
#    utm29.ImportFromEPSG(32629)  # utm 29
    utm29.ImportFromEPSG(900913)  # web mercator 3857
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(4326)
    tx = osr.CoordinateTransformation(wgs84,utm29)
    # Get the Geotransform vector
    # Work out the boundaries of the new dataset in the target projection
    (ulx29, uly29, ulz29) = tx.TransformPoint(ulx, uly)  # corner coords in utm meters
    (lrx29, lry29, lrz29) = tx.TransformPoint(ulx + xres*xsize, uly - yres*ysize )
    filenameutm = filename[0:-4] + '_web.tif'
    dest = driver.Create(filenameutm, xsize, ysize, 2, gdal.GDT_Float32)
    xres29 = round((lrx29 - ulx29)/xsize, 2) # pixel size, utm meters
    yres29 = abs(round((lry29 - uly29)/ysize, 2))
    new_gt = [ulx29, xres29, 0, uly29, 0, -yres29]
    dest.SetGeoTransform(new_gt)
    dest.SetProjection(utm29.ExportToWkt())
    gdal.ReprojectImage(ds, dest, wgs84.ExportToWkt(), utm29.ExportToWkt(), gdal.GRA_Bilinear)
    dest.GetRasterBand(1).SetNoDataValue(0.0)
    dest.GetRasterBand(2).SetNoDataValue(0.0)
    dest = None  # Flush the dataset to the disk
    ds = None  # only after the reprojected!
    print 'Image Created'

```