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, and transforming coordinates using advices from here (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: Then I used EPSG(3857) because I've read this and this 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).
With EPSG(900913) the output is georeferenced, but shifted about 3 raster cells to the north:
When I reproject the raster using ArcGIS (export in WGS_1984_Web_Mercator_Auxiliary_Sphere) the result is nearly fine:
And when I use the old code 102113 (41001,54004) the result is perfect:
The summary of my tests using all codes:
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?
yres = round(lons-lons, 4) # pixel size, degrees xres = round(lats-lats, 4) ysize = len(lats)-1 # number of pixels xsize = len(lons)-1 ulx = round(lons, 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'