17

I'm trying to perform a rather simple task.. but can't find a simple solution.. I'm trying to re-project a raster layer using GDAL in python. I've read this solution , and this , and this, and others. It seems that all of them assume some parameters like pixel size or xy origins but I need to get all of the parameters from the raster itself, it suppose to be as generic as possible, the only constant thing is the output coordinate system which is wgs84 (EPSG=4326).

I've also saw different tools to perform this task like gdalwrap or getransform, so I'm a bit confused.

To sum it all up, I need to re-project the input_raster from whatever coordinate system it is on to output_raster that will be in wgs84 coordinate system.

input_raster = r"C:\Users\B3.tif" 

output_raster = r"C:\Users\B3_proj.tif"

Could you help me in writing a python code using gdal to perform this task?

Community
  • 1
user88484
  • 1,779
  • 4
  • 26
  • 39
  • Instead of reinventing the wheel how about subprocess.popen('GDAL_Warp ... https://docs.python.org/2/library/subprocess.html http://www.gdal.org/gdalwarp.html after you've obtained the input coordinate system http://www.gdal.org/classGDALDataset.html#aa42537e1062ce254d124b29ff3ebe857 assuming it's defined, if the input raster spatial reference isn't defined you'll have to do some tricky assumptions based on your knowledge of your data. Note with GDAL_Warp if the existing spatial reference is defined one need only supply the -t_srs EPSG:4326 – Michael Stimson Mar 26 '17 at 21:45
  • 4
    Or just use gdal.Warp(output, input, dstSRS='EPSG:4326') if you have gdal >= 2.0 – user2856 Mar 26 '17 at 21:52
  • 2
    @MichaelMiles-Stimson the command is gdalwarp, not GDAL_Warp. – user2856 Mar 26 '17 at 21:53
  • Oh it's sooo monday morning @Luke! You are correct. gdal.warp( sounds like an answer to me... I use subprocess.popen because I want to do more than one concurrently usually. – Michael Stimson Mar 26 '17 at 21:55

1 Answers1

36

So, thanks to @Luke I used the simple line: gdal.Warp(output_raster,input_raster,dstSRS='EPSG:4326') and it works, this was exactly what I was looking for, simple code with few line to execute simple task.. I don't understand why all of the other answers I found online are much more complicated.. The code I wrote is:

from osgeo import gdal

filename = r"C:\path\to\input\raster
input_raster = gdal.Open(filename)
output_raster = r"C:\path\to\output\raster
warp = gdal.Warp(output_raster,input_raster,dstSRS='EPSG:4326')
warp = None # Closes the files
Colonel_Old
  • 103
  • 3
user88484
  • 1,779
  • 4
  • 26
  • 39