Raster calculator shifts raster cells during calculation

[WadeWall]

Hi all,

I was not aware of the BAI function ins Spatial Analyst, and so I had written the burned area index into a function. However, the results shifted the raster cells and I am not sure why.

Here is the example of the code I used in raster calculator:

1/((0.1 - Power("LS5_c1_19920425.tif_B3",2)) + Power("LS5_c1_19920425.tif_B4",2))

where files "..B3.tif" and "..B4.tif" are the red and NIR bands from landsat imagery.

When I use the BAI function in spatial analyst, the alignment is correct. I don't see what in the code that I have would lead to a shift in the raster cells. Any help would be appreciated.

[Luke_Pinner]

Set your snap raster and extent environments in the raster calculator.

[WadeWall]

Thanks for the suggestion, but that didn't do the trick. The issue is that the y cell size is slightly larger on the output raster, when compared to the input raster. Here is the code that I have used to try to remedy the situation:

### get and set extent of raster
file = "image.tif"

import arcpy
from arcpy import env
from arcpy import Raster
arcpy.CheckOutExtension("Spatial")
from arcpy.sa import *

arcpy.ResetEnvironments()
arcpy.env.overwriteOutput = True
arcpy.env.snapRaster  = arcpy.sa.Raster(file)
arcpy.env.cellSize = file

rast = arcpy.sa.Raster(file)
desc = arcpy.Describe(rast)
xmin = desc.extent.XMin
xmax = desc.extent.XMax
ymin = desc.extent.YMin
ymax = desc.extent.YMax

print("xmin: %s \nxmax: %s \nymin: %s \nymax: %s" % (xmin, xmax, ymin, ymax))

arcpy.env.extent = arcpy.Extent(xmin, ymin, xmax, ymax)

red = arcpy.sa.Raster(file + '\B3')
NIR = arcpy.sa.Raster(file + '\B4')
tmp1 = 0.1 - red
tmp2 = 0.06 - NIR
denom = arcpy.sa.Power(tmp1,2) + arcpy.sa.Power(tmp2,2)
BAI = arcpy.sa.Divide(1, denom)

desc = arcpy.Describe(BAI)
xmin = desc.extent.XMin
xmax = desc.extent.XMax
ymin = desc.extent.YMin
ymax = desc.extent.YMax
print("xmin: %s \nxmax: %s \nymin: %s \nymax: %s" % (xmin, xmax, ymin, ymax))

BAI.save(r"D:\test.tif")

### get cell size
cllsize = arcpy.management.GetRasterProperties(template,"CELLSIZEX")
print(cllsize.getOutput(0))

cllsize = arcpy.management.GetRasterProperties(BAI,"CELLSIZEX")
print(cllsize.getOutput(0))

cllsize = arcpy.management.GetRasterProperties(template,"CELLSIZEY")
print(cllsize.getOutput(0))

cllsize = arcpy.management.GetRasterProperties(BAI,"CELLSIZEY")
print(cllsize.getOutput(0))

Even after attempting to set the cell size and extent, the output raster is not the same as the input raster. Any other suggestions?

[WadeWall]

Okay, here is some minimal code, along with the file that is causing me problems.

import arcpy
from arcpy import env
from arcpy import Raster
arcpy.CheckOutExtension("Spatial")
from arcpy.sa import *

arcpy.ResetEnvironments()
arcpy.env.overwriteOutput = True
arcpy.env.snapRaster  = arcpy.sa.Raster(file)
arcpy.env.cellSize = file

rast = arcpy.sa.Raster(file)
desc = arcpy.Describe(rast)
xmin = desc.extent.XMin
xmax = desc.extent.XMax
ymin = desc.extent.YMin
ymax = desc.extent.YMax

print("xmin: %s \nxmax: %s \nymin: %s \nymax: %s" % (xmin, xmax, ymin, ymax))

arcpy.env.extent = arcpy.Extent(xmin, ymin, xmax, ymax)

red = arcpy.sa.Raster(file + '\B3')
NIR = arcpy.sa.Raster(file + '\B4')
tmp1 = 0.1 - red ## This is where some problems are happening
tmp2 = 0.06 - NIR

### get cell size
cllsizeX = arcpy.management.GetRasterProperties(Raster(file),"CELLSIZEX")
print("cell size X = " + cllsizeX.getOutput(0))
cllsizeY = arcpy.management.GetRasterProperties(Raster(file),"CELLSIZEY")
print("cell size Y = " + cllsizeY.getOutput(0))

### Now this is smaller than the image cell size
cllsizeX = arcpy.management.GetRasterProperties(NIR,"CELLSIZEX")
print("cell size X = " + cllsizeX.getOutput(0))
cllsizeY = arcpy.management.GetRasterProperties(NIR,"CELLSIZEY")
print("cell size Y = " + cllsizeY.getOutput(0))

### Same as NIR
cllsizeX = arcpy.management.GetRasterProperties(red,"CELLSIZEX")
print("cell size X = " + cllsizeX.getOutput(0))
cllsizeY = arcpy.management.GetRasterProperties(red,"CELLSIZEY")
print("cell size Y = " + cllsizeY.getOutput(0))

### Now this larger than NIR. GRRRR
cllsizeX = arcpy.management.GetRasterProperties(tmp1,"CELLSIZEX")
print("cell size X = " + cllsizeX.getOutput(0))
cllsizeY = arcpy.management.GetRasterProperties(tmp1,"CELLSIZEY")
print("cell size Y = " + cllsizeY.getOutput(0))

 I don't understand how the multiband raster can have cell sizes that are different than the underlying bands. This seems to be the case, but is not the case in R or QGIS. Any help would be appreciated.

[WadeWall · Accepted solution]

The issue was that I had previously used arcpy.management.Clip to trim a raster with a shapefile, and this resulted in the cell size being, not 30 m, but  something like 29.96 and 29.954. Downstream tools used different combination of these values to assign cell size and that is where the problem occurred.