Hello,
I have two rasters that were calculated using the Trend tool in the Spatial Analyst toolbox, each from a separate feature layers of points. Both rasters are linear trend surfaces with polynomial degree 1. One raster has a larger footprint than the other, and the smaller raster is completely contained within the larger raster. My goal is to extend smaller trend surface to match the footprint of the larger one so that I can subtract the two surfaces and find the difference.
According to this article the text file that is produced as part of the Trend output contains the coefficients multiplied by latitude and longitude values. I am trying to write some arcpy code (shamelessly stolen from this stack exchange answer) that will use these coefficients with the footprint and cell size of the larger raster. I would expect the output of this code to have the same values as the smaller trend surface anywhere in the smaller footprint, with the trend smoothly extended out to the larger footprint. I do get a raster with the expected footprint but the values are off by ~800 or so.
import arcpy
from arcpy.sa import *
from arcpy import env as E
arcpy.env.extent = "large_footprint"
arcpy.env.cellSize = 0.03
arcpy.gp.SingleOutputMapAlgebra_sa("-138.163214365158 +
0.00280100650986114 * $$ROWMAP +
-0.00594606008855511 * $$COLMAP",
"extended_trend")
I am not married to this method - if anyone has a different suggestion I would love to try it. I am not sure how to get around the scaling up of the footprint without somehow referencing the columns and rows of the larger footprint.
I am attaching the text file produced by Trend that has the coefficients listed above. If there is a way I can attach .tif files of the original and extended trend surfaces please let me know. In the meantime I also posted a screenshot of the original and extended surfaces in the contents panel to show the differences (~600 for original trend, ~-140 for the extended surface).
Solved! Go to Solution.
I solved the issue. My error came from my misunderstanding of $$ROWMAP and $$COLMAP. These do not directly return the easting and northing of the raster cells, but an incremented list of the rows and columns of the raster, i.e. [1, 2, ..., n] for both rows and columns. (In hindsight I probably should have realized this sooner...)
Below is the code I used to generate the extended trend surface I needed.
import numpy as np
import arcpy
from arcpy.sa import *
from arcpy import env as E
arcpy.env.extent = "large_footprint"
arcpy.env.cellSize = 0.3
cellSize = float(arcpy.env.cellSize)
xmin = arcpy.env.extent.XMin
xmax = arcpy.env.extent.XMax
ymin = arcpy.env.extent.YMin
ymax = arcpy.env.extent.YMax
easting = np.arange(xmin, xmax, cellSize)[:-1] #trim due to double representation
northing = np.arange(ymin, ymax, cellSize)[::-1] #need reverse order
xx, yy = np.meshgrid(easting, northing)
trend = -138.163214365158 + 0.00280100650986114 * xx + -0.00594606008855511 * yy
trend_extended = arcpy.NumPyArrayToRaster(trend, arcpy.Point(xmin, ymin), x_cell_size=cellSize)
trend_extended.save("trend_extended")
I solved the issue. My error came from my misunderstanding of $$ROWMAP and $$COLMAP. These do not directly return the easting and northing of the raster cells, but an incremented list of the rows and columns of the raster, i.e. [1, 2, ..., n] for both rows and columns. (In hindsight I probably should have realized this sooner...)
Below is the code I used to generate the extended trend surface I needed.
import numpy as np
import arcpy
from arcpy.sa import *
from arcpy import env as E
arcpy.env.extent = "large_footprint"
arcpy.env.cellSize = 0.3
cellSize = float(arcpy.env.cellSize)
xmin = arcpy.env.extent.XMin
xmax = arcpy.env.extent.XMax
ymin = arcpy.env.extent.YMin
ymax = arcpy.env.extent.YMax
easting = np.arange(xmin, xmax, cellSize)[:-1] #trim due to double representation
northing = np.arange(ymin, ymax, cellSize)[::-1] #need reverse order
xx, yy = np.meshgrid(easting, northing)
trend = -138.163214365158 + 0.00280100650986114 * xx + -0.00594606008855511 * yy
trend_extended = arcpy.NumPyArrayToRaster(trend, arcpy.Point(xmin, ymin), x_cell_size=cellSize)
trend_extended.save("trend_extended")