2

The equivalent question has been asked wrt ArcMap, but I'm using QGIS and am also an inexperienced user.

I have multiple rasters, provided as GeoTIFFs, that are all projected to the same custom CRS but are often offset by a few meters. I would like to treat one of them as correct and adjust the other rasters to match visually. A translation (shift by x/y) is all that's needed.

Is there a simple WYSIWYG-approach using QGIS to manually/visually translate one raster relative to its CRS or relative to another raster (and its CRS)?

It would be possible to georeference the rasters again, but they are quite close already, and I'm concerned about introducing errors of rotation and scale by doing that. I also don't know yet how to do that. I'm searching first for an answer that starts from the current, very close georeferencing solution.

If a WYSIWYG-approach doesn't exist, I could repeatedly guess at the translations necessary if I could modify the coordinates listed under Extent on the Information page of the raster's property sheet. Is that possible in QGIS?

Failing that, is there a way to translate the images using GDAL command line utilities or Python using GDAL?

The answer for the ArcGIS version of this question mentioned a solution using a datum shift and another approach for raster data using ArcGIS' georeferencing toolbar. I haven't found equivalent functionality in QGIS yet.

Improve this question
2

1 Answer 1

Reset to default
4

There may certainly be an easier/ better method, but you can try this script using the GDAL Python API, based on some of the examples found here.

You will need to edit first four executed lines after the import statement to define the paths to your input raster, path to where you want the output raster saved, then define your X and Y offset values (which may be positive of negative and are in the input raster CRS units).

You can run this script from an editor in the QGIS Python console.

if you have multiple rasters, all with the same required offsets, it would be easy enough to modify the script to loop over all .tif files in a directory.

from osgeo import gdal, osr

# Define input and output file paths
in_path = 'C:\\Path\\To\\Input_raster.tif'
out_path = 'C:\\Path\\To\\Output_raster.tif' # Will be created!

# Define desired X and Y offsets in input raster Spatial Reference System (SRS) units...
# ...using an accurate meter-based SRS e.g. an appropriate UTM CRS is strongly recommended
# Can be + or - values
x_offset = 250
y_offset = 250

raster = gdal.Open(in_path)

# Get input raster bands
input_bands = [raster.GetRasterBand(i+1) for i in range(raster.RasterCount)]

# Get input raster data type
data_type = input_bands[0].DataType

# Create list of bands read into arrays
band_arrays = [band.ReadAsArray() for band in input_bands]

geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize

driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(out_path, cols, rows, len(input_bands), data_type)
outRaster.SetGeoTransform((originX+x_offset, pixelWidth, 0, originY+y_offset, 0, pixelHeight))

# Write array from each input band to each output band
for index, array in enumerate(band_arrays):
    outband = outRaster.GetRasterBand(index+1)
    if outband:
        outband.WriteArray(array)
        outband.FlushCache()
    
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
outRaster.SetProjection(outRasterSRS.ExportToWkt())

outRaster = None

Result on a sample raster dataset:

enter image description here

Edit:

I have written a script using the QgsMapToolEmitPoint class which may improve your workflow if you can visually locate the same location/pixel on two slightly misaligned rasters. When you run this tool you can left-click on the raster you wish to move, then left-click again at the same location on your reference raster. The two points will be marked by QgsRubberBand point markers. You can re-click the two reference points until you are happy with them. Then right-click to finish. The two points will be used to calculate the X and Y offset values.

Paste the code below into an editor in the QGIS Python console. Select the layer you wish to move as the active layer, and click run.

Full script:

from osgeo import gdal, osr

class MapToolShiftRaster(QgsMapToolEmitPoint):
    
    def __init__(self, canvas, shift_layer):
        self.canvas = canvas
        self.shift_layer = shift_layer
        QgsMapToolEmitPoint.__init__(self, self.canvas)
        self.first_click = True
        self.pnt_A = None
        self.pnt_B = None
        self.X_offset = None
        self.Y_offset = None
        
        self.rb1 = None
        self.rb2 = None
        self.rb3 = None
    
    def canvasMoveEvent(self, e):
        if not self.first_click:
            if self.rb2:
                self.rb2.reset()
            geom = self.transform_geom(QgsGeometry().fromPolylineXY([self.pnt_A, self.toLayerCoordinates(self.shift_layer, e.mapPoint())]))
            self.rb2 = QgsRubberBand(self.canvas, QgsWkbTypes.LineGeometry)
            self.rb2.setToGeometry(geom)
            self.rb2.setColor(QColor('Blue'))
            self.rb2.show()
    
    def canvasReleaseEvent(self, e):
        if e.button() == Qt.LeftButton:
            if self.first_click:
                if self.rb1:
                    self.rb1.reset()
                if self.rb2:
                    self.rb2.reset()
                if self.rb3:
                    self.rb3.reset()
                self.pnt_A = self.toLayerCoordinates(self.shift_layer, e.mapPoint())
                geom = self.transform_geom(QgsGeometry().fromPointXY(self.pnt_A))
                self.rb1 = QgsRubberBand(self.canvas, QgsWkbTypes.PointGeometry)
                self.rb1.setToGeometry(geom)
                self.rb1.setIcon(QgsRubberBand.ICON_BOX)
                self.rb1.setColor(QColor('Red'))
                self.rb1.show()
                self.first_click = False
            elif not self.first_click:
                if self.rb2:
                    self.rb2.reset()
                self.pnt_B = self.toLayerCoordinates(self.shift_layer, e.mapPoint())
                if self.pnt_A and self.pnt_B:
                    geom = self.transform_geom(QgsGeometry().fromPointXY(self.pnt_B))
                    self.rb3 = QgsRubberBand(self.canvas, QgsWkbTypes.PointGeometry)
                    self.rb3.setToGeometry(geom)
                    self.rb3.setIcon(QgsRubberBand.ICON_BOX)
                    self.rb3.setColor(QColor('Red'))
                    self.rb3.show()
                    self.X_offset = self.pnt_B.x() - self.pnt_A.x()
                    self.Y_offset = self.pnt_B.y() - self.pnt_A.y()
                self.first_click = True
        elif e.button() == Qt.RightButton:
            if self.rb1:
                self.rb1.reset()
            if self.rb2:
                self.rb2.reset()
            if self.rb3:
                self.rb3.reset()
            self.shift_raster(self.shift_layer, self.X_offset, self.Y_offset)
            
#########################################################################

    def shift_raster(self, lyr, x_offset, y_offset):
        if x_offset is not None and y_offset is not None:
            in_path = lyr.source()
            print(in_path)
            out_string = QFileDialog.getSaveFileName(None, 'Save output raster', filter='.tif')
            file_name = out_string[0].split('/')[-1]
            out_path = ''.join(out_string)

            raster = gdal.Open(in_path)
            if not raster:
                print('Input raster is none')
                return

            # Get input raster bands
            input_bands = [raster.GetRasterBand(i+1) for i in range(raster.RasterCount)]

            # Get input raster data type
            data_type = input_bands[0].DataType

            # Create list of bands read into arrays
            band_arrays = [band.ReadAsArray() for band in input_bands]

            geotransform = raster.GetGeoTransform()
            originX = geotransform[0]
            originY = geotransform[3]
            pixelWidth = geotransform[1]
            pixelHeight = geotransform[5]
            cols = raster.RasterXSize
            rows = raster.RasterYSize

            driver = gdal.GetDriverByName('GTiff')
            outRaster = driver.Create(out_path, cols, rows, len(input_bands), data_type)
            if not outRaster:
                print('Output raster is none')
                return
            outRaster.SetGeoTransform((originX+x_offset, pixelWidth, 0, originY+y_offset, 0, pixelHeight))

            # Write array from each input band to each output band
            for index, array in enumerate(band_arrays):
                outband = outRaster.GetRasterBand(index+1)
                if outband:
                    outband.WriteArray(array)
                    outband.FlushCache()
                
            outRasterSRS = osr.SpatialReference()
            outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
            outRaster.SetProjection(outRasterSRS.ExportToWkt())

            outRaster = None
            
            iface.addRasterLayer(out_path, file_name, 'gdal')

#########################################################################
            
    def transform_geom(self, geom):
        proj_crs = QgsProject.instance().crs()
        lyr_crs = self.shift_layer.crs()
        if proj_crs != lyr_crs:
            xform = QgsCoordinateTransform(lyr_crs, proj_crs, QgsProject.instance())
            geom.transform(xform)
            return geom
        return geom
            
    def deactivate(self):
        if self.rb1:
            self.rb1.reset()
            self.rb1 = None
        if self.rb2:
            self.rb2.reset()
            self.rb2 = None
        if self.rb3:
            self.rb3.reset()
            self.rb3 = None
        

layer = iface.activeLayer()
if layer.type() == QgsMapLayerType.RasterLayer:
    T = MapToolShiftRaster(iface.mapCanvas(), layer)
    iface.mapCanvas().setMapTool(T)
else:
    iface.messageBar().pushMessage('Please select a raster layer!')

See the screencast below for a short demo:

enter image description here

You can see the results after one run below. I was able to reduce the misalignment from over 250m to around 0.6m. If your rasters are only a few meters apart you should be able to align them almost perfectly in one go by zooming right in. Otherwise you should be able to achieve a satisfactory result on the 2nd go.

enter image description here

Improve this answer
2
  • Thank you very much. This does work to shift the rasters. I ended up writing almost the same code with the only difference being a call to shutil.copyfile() first followed by changing the geotransform on the new file in place. Your method is better, and I'll change my code. To use this solution, I repeatedly loaded and unloaded my rasters, guessing at the offsets and re-running the python code. I'm still hoping for a solution within QGIS that involves dragging the raster in some way that gives immediate feedback as the user matches features visually. That would be much quicker.
    – mshirley
    Commented Nov 9, 2021 at 17:24
  • @mshirley, I have updated my answer with a script which uses the QgsMapToolEmitPoint subclass to give you a somewhat visual approach to aligning rasters. I hope it helps you. Let me know how it goes.
    – Ben W
    Commented Nov 10, 2021 at 10:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.