ArcGIS / Python / Práce s rastrovými daty: Porovnání verzí
Řádek 32: | Řádek 32: | ||
== Fokální operace -moving window (kernel) == | == Fokální operace -moving window (kernel) == | ||
<source lang=python> | |||
import arcpy | |||
import numpy | |||
from scipy.ndimage import filters | |||
# promenne prostredi | |||
arcpy.env.workspace = "T:\\Groups\\FPF\\spatial\\IDU\\Deschutes2" | |||
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 == | == Odkazy == |
Verze z 19. 12. 2013, 16:11
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 cvičení předmětu ZODH.
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")
Fokální operace -moving window (kernel)
import arcpy
import numpy
from scipy.ndimage import filters
# promenne prostredi
arcpy.env.workspace = "T:\\Groups\\FPF\\spatial\\IDU\\Deschutes2"
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")