Zobrazení zdroje stránky ArcGIS / Python / Práce s rastrovými daty

{{upravit}}
{{ArcGIS}}

== Výpočet normalizovaného diferenčního vegetačního indexu ==

Výpočet normalizovaného diferenčního vegetačního indexu (NDVI) - viz [[153ZODH / 4. cvičení#NDVI|cvičení předmětu ZODH]].

<source lang=python>
import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")

arcpy.env.overwriteOutput = True
arcpy.env.workspace = "c:/Users/landa/arcpy" 

### tm3_name = "tm3.tif"
### tm4_name = "tm4.tif"
tm3_name = arcpy.GetParameterAsText(0)
tm4_name = arcpy.GetParameterAsText(1)

# nacteni vstupnich dat
tm3 = Raster(tm3_name)
tm4 = Raster(tm4_name)

# vypocet NDVI
ndvi = 1.0 * (tm4 - tm3) / (tm4 + tm3)

# ulozeni vysledku do souboru TIFF
ndvi.save("ndvi.tif")
</source>

== Fokální operace - moving window (kernel) ==

Požadavky:

* [http://www.numpy.org/ NumPy]
* [http://scipy.org/scipylib/index.html SciPy library]

<source lang=python>
import arcpy
import numpy
from scipy.ndimage import filters

# promenne prostredi
arcpy.env.workspace = "c:/Users/landa/arcpy"
rasterCellSize = 30
rasterInputName = "tm3.tif"
rasterOutputName = 'tm3_filtered'
arcpy.env.cellSize = rasterCellSize

kernelFP = numpy.array([[False,True,False], [True,True,True], [False,True,False]])

# funkce kernelu
def neighborFunc(kernelArray):
    if (kernelArray[2] != -999):
        return kernelArray[0] + kernelArray[1] + kernelArray[3] + kernelArray[4]
    else:
        return kernelArray[2]

# popis vstupu
rasterDesc = arcpy.Describe(rasterInputName)
lowerLeftCorner = arcpy.Point(rasterDesc.Extent.XMin, rasterDesc.Extent.YMin)

# nahrani vstupni rastrove vrstvy do pole numpy
myArray = arcpy.RasterToNumPyArray(rasterInputName, lowerLeftCorner, rasterDesc.width, rasterDesc.height, -999)

# aplikace moving window
newRaster = filters.generic_filter(myArray, neighborFunc, footprint=kernelFP)

# ulozeni vysledku do nove rastrove vrstvy
rasterToSave = arcpy.NumPyArrayToRaster(newRaster, lowerLeftCorner, rasterCellSize, rasterCellSize, -999)
rasterToSave.save(rasterOutputName)
arcpy.DefineProjection_management(rasterOutputName, rasterDesc.spatialReference)
arcpy.BuildRasterAttributeTable_management(rasterOutputName, "Overwrite")
</source>

== Odkazy ==

* http://blog.remotesensing.io/2013/05/using-arcpy-for-raster-analysis/
* http://sdmg.forestry.oregonstate.edu/olsen-ex1