21

I am VERY new to this whole GIS thing.

I am trying to convert some projected geoTiff images to WGS84 using gdal in python. I have found a post that outlines the process to transform points within projected GeoTiffs using something similar to the following:

from osgeo import osr, gdal

get the existing coordinate system


ds = gdal.Open('path/to/file')
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())


create the new coordinate system


wgs84_wkt = """
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]"""
new_cs = osr.SpatialReference()
new_cs .ImportFromWkt(wgs84_wkt)


create a transform object to convert between coordinate systems


transform = osr.CoordinateTransformation(old_cs,new_cs)


#get the point to transform, pixel (0,0) in this case
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
miny = gt[3] + widthgt[4] + heightgt[5]


#get the coordinates in lat long
latlong = transform.TransformPoint(x,y)

If I want to convert these points and create a new WGS84 GeoTiff file, is this the best way to go about it?

Is there a function that exists that will do such as task in 1 step?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
JT.WK
  • 313
  • 1
  • 2
  • 5

2 Answers2

26

One simpler way would be to use the GDAL command line tools:

gdalwarp infile.tif outfile.tif -t_srs "+proj=longlat +ellps=WGS84"

That can be invoked easily enough via scripting for batch jobs. This will warp the raster to a new grid, and there are other options to control the details.

Link

The target (t_srs) coordinate system can be PROJ.4 as shown, or an actual file with WKT, amongst other options. The source (-s_srs) coordinate system is assumed known, but can be set explicitly like the target.

That might be easier to get the job done if you don't have to use the GDAL API via Python.

There is a tutorial for warping with the API here, it says the support for Python is incomplete (I don't know the details): Link

Glorfindel
  • 1,096
  • 2
  • 9
  • 14
mdsumner
  • 8,196
  • 1
  • 27
  • 30
17

As mdsumner said, it's much easier to use command line than the python bindings, unless you want to execute very complex tasks.

So, if you like python, like I do, you can run the command line tool with:

import os  
os.sys('gdalwarp infile.tif outfile.tif -t_srs "+proj=longlat +ellps=WGS84"')

or iterate through a list of files:

listoffiles=['file1, file2, file3, filen']  
for file in listoffiles:
    os.system('gdalwarp %s %s -t_srs "+proj=longlat +ellps=WGS84"' %(file,'new_'+file))

And even use multiprocessing tools use the full power of your machine to execute big tasks.

nickves
  • 11,519
  • 3
  • 42
  • 75
Pablo
  • 9,827
  • 6
  • 43
  • 73