2

I want to crop a landcover raster (cell size 10 meters) to my project region. The problem is that the new raster has almost always an offset.

In the example I show below I made the new raster transparent: one can see the offset, which is more than 4 meters.

In ArcGIS there is an option in the toolbox environment to snap to another raster. Does it exist in QGIS or what could be the workaround?

I found another post here Specifying snap raster in QGIS? but its 2,5 years old, maybe it is possible nowadays. I know it depends on the coordiantes for the extent. But finding the right extent coordinates in order to prevent the offset is quite time consuming.

enter image description here

Enzo Baldini
  • 2,327
  • 1
  • 15
  • 29

1 Answers1

1

I wrote some code, which is the best solution I'm aware of at this stage when nice alternatives aren't available (e.g having pixel coordinates that are all divisible by 10cm or something)

import numpy
from qgis.core import QgsRasterLayer
import os
import glob

"""
##############################################################
User options
"""


#The raster that is to be snapped
inputRaster     = 'D:/Imagery/ToSnap.tif'


#The raster to snap to
snapRaster      = 'D:/Imagery/OriginalBaseRaster.tif'


#The output raster to be created
outputRaster    = 'D:/Imagery/Snapped.tif'


#Options for output compression
compressOptions = 'COMPRESS=LZW|PREDICTOR=1|NUM_THREADS=ALL_CPUS|BIGTIFF=IF_SAFER|TILED=YES'


"""
##############################################################
Input raster info
"""


inputRas = QgsRasterLayer(inputRaster) 
pixelSizeXInputRas = inputRas.rasterUnitsPerPixelX()
pixelSizeYInputRas = inputRas.rasterUnitsPerPixelY()


inputRasBounds = inputRas.extent() 
xminInputRas = inputRasBounds.xMinimum()
xmaxInputRas = inputRasBounds.xMaximum()
yminInputRas = inputRasBounds.yMinimum()
ymaxInputRas = inputRasBounds.yMaximum()
coordsInputRas = "%f %f %f %f" %(xminInputRas, ymaxInputRas, xmaxInputRas, yminInputRas)
print(coordsInputRas)


"""
##############################################################
Snap raster info
"""


snapRas = QgsRasterLayer(snapRaster) 
pixelSizeXSnapRas = snapRas.rasterUnitsPerPixelX()
pixelSizeYSnapRas = snapRas.rasterUnitsPerPixelY()
coordinateSystemSnapRas = snapRas.crs().authid()
snapRasterHeightSnapRas = snapRas.height()
snapRasterWidthSnapRas = snapRas.width()


snapRasBounds = snapRas.extent()
xminSnapRas = snapRasBounds.xMinimum()
xmaxSnapRas = snapRasBounds.xMaximum()
yminSnapRas = snapRasBounds.yMinimum()
ymaxSnapRas = snapRasBounds.yMaximum()
coordsSnapRas = "%f %f %f %f" %(xminSnapRas, ymaxSnapRas, xmaxSnapRas, yminSnapRas)
print(coordsSnapRas)


"""
##############################################################
Get a list of possible snap coordinates
"""


listOfXCoords = []
nextXCoord = xminSnapRas


while nextXCoord < xmaxSnapRas:
    nextXCoord = nextXCoord + pixelSizeXSnapRas
    listOfXCoords.append(nextXCoord)
print(str(listOfXCoords[:10]) + '...')


listOfYCoords = []
nextYCoord = yminSnapRas


while nextYCoord < ymaxSnapRas:
    nextYCoord = nextYCoord + pixelSizeYSnapRas
    listOfYCoords.append(nextYCoord)
print(str(listOfYCoords[:10]) + '...')


"""
##############################################################
Determine what the delta X & Y needs to be
"""


#We'll aim to snap the top left corner (highest y val, lowest x val)
closestXCoordinate = min(listOfXCoords, key=lambda x:abs(x-xminInputRas))
closestYCoordinate = min(listOfYCoords, key=lambda x:abs(x-ymaxInputRas))


shiftX = closestXCoordinate - xminInputRas
shiftY = closestYCoordinate - ymaxInputRas


shiftParameter = '-a_ullr ' + str(xminInputRas + shiftX) + ' ' + str(ymaxInputRas + shiftY) + ' ' + str(xmaxInputRas + shiftX) + ' ' + str(yminInputRas + shiftY)


processing.run("gdal:translate", {'INPUT':inputRaster,'TARGET_CRS':None,'NODATA':None,'COPY_SUBDATASETS':False,'OPTIONS':compressOptions,'EXTRA':shiftParameter,'DATA_TYPE':0,'OUTPUT':outputRaster})

This determines the closest pixel corner in a base raster to the top left corner of an input raster, then puts the shifted extent in the parameter -a_ullr in gdal_translate

This won't work when the top left corner of the input raster is well beyond the extent of the base raster. It also requires the two rasters to have the same CRS. You could modify the code to account for this if you wanted I guess

If there are any improvements to this let me know

Buff Fox
  • 119
  • 11