1 / 11

Raster MANIPULATION IN PYTHON

Raster MANIPULATION IN PYTHON. Ditch ArcGIS and get to know NUMPY and SCIPY. What we need to start. ArcGIS 10x Python 2.7.2 or better Numpy (installed by ArcGIS ) SciPy Some luck. RasterToNumpyArray. Check projection consistency Find a corner point based on a reference raster

tarala
Télécharger la présentation

Raster MANIPULATION IN PYTHON

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Raster MANIPULATION IN PYTHON Ditch ArcGIS and get to know NUMPY and SCIPY

  2. What we need to start • ArcGIS 10x • Python 2.7.2 or better • Numpy (installed by ArcGIS) • SciPy • Some luck

  3. RasterToNumpyArray • Check projection consistency • Find a corner point based on a reference raster • Determine number of rows and columns to load • Load raster data with RasterToNumpyArray

  4. Masked Arrays • Write raster image as masked array. • Why masked array? • >>> import numpy.ma as ma • >>> x = [0.,1.,-9999.,3.,4.] • >>> mx = ma.masked_values (x, -9999.) • >>> print mx.mean() • 2.0

  5. Custom Filters with scipy • What is scipy? SciPy (pronounced “Sigh Pie”) is a Python-based ecosystem of open-source software for mathematics, science, and engineering. • Ndimage. Multidimentional image processing. http://docs.scipy.org/doc/scipy/reference/ndimage.html • How does one build a custom filter? Generic_filter

  6. Scipy – generic_filter • Input – Input array to filter • Function – Function to apply at each element • Footprint - footprint is a boolean array that specifies (implicitly) a shape, but also which of the elements within this shape will get passed to the filter function. • Mode - {‘reflect’, ‘constant’, ‘nearest’, ‘mirror’, ‘wrap’} • Cval - Value to fill past edges of input if mode is ‘constant’. Default is 0.0

  7. Boundary Conditions • Since the output elements are a function of elements in the neighborhood of the input elements, the borders of the array need to be dealt with appropriately by providing the values outside the borders. This is done by assuming that the arrays are extended beyond their boundaries according certain boundary conditions. In the functions described below, the boundary conditions can be selected using the mode parameter which must be a string with the name of the boundary condition. Following boundary conditions are currently supported: • “nearest” Use the value at the boundary [1 2 3]->[1 1 2 3 3] • “wrap” Periodically replicate the array [1 2 3]->[3 1 2 3 1] • “reflect” Reflect the array at the boundary [1 2 3]->[1 1 2 3 3] • “constant” Use a constant value, default is 0.0 [1 2 3]->[0 1 2 3 0] • The “constant” mode is special since it needs an additional parameter to specify the constant value that should be used.

  8. How to build a footprint • def buildCircularKernel(radiusMeters, referenceDesc, logFile): • kernel = np.array([], dtype=bool) • # build kernel from radius and input raster cell size • window = int(radiusMeters / referenceDesc.meanCellHeight) • for row in range(-window, window + 1): • for col in range(-window, window + 1): • Xdist = abs(col * referenceDesc.meanCellHeight) • Ydist = abs(row * referenceDesc.meanCellHeight) • totalDist = pow(Xdist**2 + Ydist**2, 0.5) • if totalDist > radiusMeters: • kernel = np.append(kernel, [False]) • else: • kernel = np.append(kernel, [True]) • kernel.shape = (pow(kernel.size,0.5),pow(kernel.size,0.5)) • print str(kernel.size) + " pixels in circular kernel" • return(kernel) • OR… • kernelFP = np.array([[False, True, False], [True, True, True], [False, True, False]], dtype=bool)

  9. Writing a Custom Function • # function to generate average value in kernel that is not NODATA • def meanFunc(kernel): • if kernel[int(kernel.size / 2)] == nodataVal: • return nodataVal • else: • # remove backgroud values from array • return np.average(kernel[np.where( kernel != nodataVal )])

  10. A more complex filter • # function to generate edge contrast index in kernel that is not NODATA • def edgeContrastFunc(kernel, edgeDict, argsFilter): • # if center pixel is bigger than 'early' threshold, skip • # nodata pixels are still scanned in case they have large neighbor which would indicate an edge • # for this exercise, nodata is assumed to be open or 'early' • if (kernel[int(kernel.size / 2)] > argsFilter[2]): • return nodataVal • else: • # search early pixel for edges and record in early pixel dictionary • edgeContrastList = [] • for adjPixel in (0, 1, 3 ,4): • if (kernel[adjPixel] > argsFilter[3]): • edgeContrast = 1.0 • elif (kernel[adjPixel] > argsFilter[2]): • edgeContrast = (argsFilter[4] * kernel[adjPixel]) + argsFilter[5] • else: • edgeContrast = -1.0 • # edge found! add to list that will be inserted into the dictionary and calc stats • if (edgeContrast > -1.0): • edgeContrastList.append(edgeContrast) • argsFilter[0] += 1 • if edgeContrast > argsFilter[1]: • argsFilter[1] = edgeContrast • if (len(edgeContrastList) == 0): • return nodataVal • else: • earlyPatchNum = len(edgeDict) + 1 • edgeDict[earlyPatchNum] = edgeContrastList • return earlyPatchNum

  11. Questions?? • http://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.filters.generic_filter.html#scipy.ndimage.filters.generic_filter

More Related