7

I'd like to use differently sized windows as I process (a.k.a. filter) a raster, with the size varying according to latitude.

There are tools that allow a user to specify custom windows (e.g., "User-Defined Filter" in SAGA), but I don't see a way to dynamically adjust the size of the window according to some other variable (e.g., the values in a second, coincident raster which could easily simply state their latitudes).

Finding local maxima with variable window search using ArcGIS Desktop? addresses this with respect to discrete areas; instead, I'd like to be able to continuously vary the size of the window.

Can anyone think of a way to do this?

Ideally in SAGA, QGIS, or GDAL. Also open to doing it in Python (e.g., numpy arrays), R (matrices), or ArcGIS. Programmatically would be ideal, but something like Raster Calculator could work too.

  • note that if your project your data in a appropriate (custom) coordinate system, a fixed window in this CRS could correspond with a variable size window in the non distorted CRS. – radouxju Mar 12 '16 at 12:33
4
+50

Input square window size raster:

enter image description here

Script:

import arcpy, os, traceback, sys, numpy
from arcpy import env
env.overwriteOutput = True

dem=r'D:\Scratch\demG'
wsize=r'D:\Scratch\wsize'
d=arcpy.Describe(dem)
try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)

    demArray = arcpy.RasterToNumPyArray(dem,"","","",-9999)
    sizeArray=arcpy.RasterToNumPyArray(wsize,"","","",-9999)
    blankArray=arcpy.RasterToNumPyArray(dem,"","","",-9999)
    d=arcpy.Describe(wsize)
    origin=d.extent.lowerLeft
    cSize=arcpy.Raster(wsize).meanCellHeight
    nRows,nCols=demArray.shape
    cellsTotal=nCols*nRows

#  main loop
    arcpy.SetProgressor("step", "", 0, nRows)
    for nRow in range (nRows):
            for nCol in range (nCols):
                S=sizeArray[nRow,nCol]
                shift=S/2
                top=max(0,nRow-shift)
                bottom=min(nRows-1,nRow+shift)
                rows=demArray[top:bottom+1]
                left=max(0,nCol-shift)
                right=min(nCols-1,nCol+shift)
                kernel=rows[:,range(left,right+1)]
                nR,nC=kernel.shape
                blankArray[nRow,nCol]=float(numpy.sum(kernel))/nR/nC
            arcpy.SetProgressorPosition()
    myRaster = arcpy.NumPyArrayToRaster(blankArray,origin,cSize,cSize)
    oneGrid=Con(myRaster<>-9999,myRaster)
    oneGrid.save(r'D:\Aerials\nztm\toriver')
    del demArray,blankArray,sizeArray
except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()        

Hillshade before:

enter image description here

and after

enter image description here

  • Nice @FelixIP, thanks :) I can see how this can be done outside of Arc too, since the essential logic is really on numpy arrays. – Paulo Raposo Mar 18 '16 at 0:24
1

Very interesting, for a fundamental approach, you could use OpenCV >= version 3.0, to combine advanced image processing routines with GDAL as a data provider.

  1. OpenCV can use GDAL as an image data provider (sorry, I found only a C++ example), if you've compiled OpenCV with GDAL support.

  2. Both libraries have python bindings.

  3. It should be easy, to import GDAL and OpenCV in the same script to write your image processing algorithm with different window sizes in OpenCV, feeded by the pixel to world calculations of GDAL.

Traditionally and for performance reasons, is the size of a convolutional image processing kernel fixed while number crunching the stuff, but a line by line operation in this context should work.

1

ArcGIS Spatial Analyst has an irregular shaped neighborhood that can be used with it's Block and Focal statistics, which can be used for filtering. Kernel files (ASCII text) can be created and implemented programmatically with Python. The linked page provides some examples.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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