Wednesday 19 October 2016

How to translate (reposition) a raster in Python?


I tried to use QGIS GeoReferencer to translate a raster layer, but somehow failed. It just goes to somewhere unexpected.


Now I decide to do it myself, with python. I think some libraries like GDAL must have done it. I only want to do an x-offset and y-offset with the tif, with no scale and rotation. Do you know what module I could use?



Answer



Depending on the raster format, you can either edit the world file, or use GDAL/Python:


from osgeo import gdal


# Open in read/write mode
rast_src = gdal.Open(rast_fname, 1)

# Get affine transform coefficients
gt = rast_src.GetGeoTransform()
# (2776450.0, 100.0, 0.0, 6352650.0, 0.0, -100.0)

The geotransform gt object is a 6-parameter tuple described here. You want to edit the first and fourth items:


# Convert tuple to list, so we can modify it
gtl = list(gt)

gtl[0] += 1000.0 # Move east 1 km
gtl[3] -= 20000.0 # Move south 20 km
# [2777450.0, 100.0, 0.0, 6332650.0, 0.0, -100.0]

# Save the geotransform to the raster
rast_src.SetGeoTransform(tuple(gtl))
rast_src = None # equivalent to save/close

Note: I tested this with QGIS 1.7. You can keep the previous raster in during this whole operation, but it stays in place the whole time. To see the updated change you can either (1) close/open the QGIS project, or (2) re-add the same raster. The second option has the advantage of being able to see the differences between "before" and "after" transform operations.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...