Home > python > Reprojecting WGS 1984 Web Mercator (EPSG:3857) in Python with GDAL

Reprojecting WGS 1984 Web Mercator (EPSG:3857) in Python with GDAL

January 21Hits:9
Advertisement

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: Reprojecting WGS 1984 Web Mercator (EPSG:3857) in Python with GDAL
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). Reprojecting WGS 1984 Web Mercator (EPSG:3857) in Python with GDAL

With EPSG(900913) the output is georeferenced, but shifted about 3 raster cells to the north: Reprojecting WGS 1984 Web Mercator (EPSG:3857) in Python with GDAL

When I reproject the raster using ArcGIS (export in WGS_1984_Web_Mercator_Auxiliary_Sphere) the result is nearly fine: Reprojecting WGS 1984 Web Mercator (EPSG:3857) in Python with GDAL

And when I use the old code 102113 (41001,54004) the result is perfect: Reprojecting WGS 1984 Web Mercator (EPSG:3857) in Python with GDAL

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?

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' 

Answers

I would reproject the files with gdalwarp.

I've done the same for files in EPSG:3763 that I want to convert to EPSG:3857. I compared the results using QGIS and Geoserver and the generated images were fine. Since a small rotation is applied to the images, you might get some black lines on the border (but these lines can be made transparent afterwards).

Since you have several tif images, you can use a script like this that does not change any existing file, and puts the generated files in a folder called 3857:

#!/bin/bash
mkdir 3857
for file in $(ls *.tif); do
    gdalwarp -s_srs EPSG:3763 -t_srs EPSG:3857 $file 3857/$file;
    listgeo -tfw 3857/$file;
done

If you also want to generate the .twf files, I've added listgeo.

This script is for Linux, but you can write something similar for Windows.

Would go for GDALwarp as well. Be sure to be consistent with "postings" and "cell" interpretations of rasters. http://www.remotesensing.org/geotiff/spec/geotiff2.5.html

Related Articles

Copyright (C) 2017 ceus-now.com, All Rights Reserved. webmaster#ceus-now.com 14 q. 0.718 s.