I am trying to make a custom tool to Arcmap using Python. It is supposed to read in a raster DEM, a point shape file and a distance, and output another raster.
The output raster should be essentially the same as the input one. In the process however, I want to edit some the raster values according to a specified function, namely those that are within the input distance from a point in the point file.
The problem with my script is the saving part. I have two different approaches. Alt 1 is working in the meaning that I get an output file. However is the projection and transformation completely off and I'm not sure how to fix it. Alt 2 is my preferred approach, but i can't figure out how to write my array to a band in the file. Also I am not sure that the edits gets saved back to the array properly, at least they are not visible in the files I do get.
Python is still a bit new to me, so there might just be something simple that I've missed. I'm grateful for all the help i can get!
I'm using Arcmap 10.0 and Python 2.6 on Windows 7.
import os, sys, json, osgeo, arcpy, math, arcgisscripting, numpy
from osgeo import gdal, osr, ogr, gdal_array
from osgeo.gdalconst import *
inDEM = arcpy.GetParameterAsText(0) # .tif
Outlets = arcpy.GetParameterAsText(1) # Point shapefile
outDEM = arcpy.GetParameterAsText(2) # Path
buffr = arcpy.GetParameterAsText(3) # Value
# Register GDAL drivers
gdal.AllRegister()
## Open the DEM
dem = gdal.Open(inDEM, GA_ReadOnly)
band1d = dem.GetRasterBand(1)
cols = dem.RasterXSize
rows = dem.RasterYSize
trf = dem.GetGeoTransform()
prj = dem.GetProjection()
driver = dem.GetDriver()
# Get info
band_info = arcpy.Describe(inDEM + '/Band_1')
nodata_dem = band_info.noDataValue
cellsz = band_info.meanCellHeight
area = int(buffr[0])
buffers = area/cellsz # The area around each chosen cell to use for calculations
# Convert to point file to raster
outlet_raster = 'testraster'
arcpy.env.overwriteOutput = True
arcpy.env.extent = inDEM
arcpy.PointToRaster_conversion(Outlets, "Enabled", outlet_raster, "MAXIMUM", "", cellsz)
arcpy.env.overwriteOutput = False
# Convert rasters to arrays (can probably be done in a more memory-efficient way)
outl_arr = arcpy.RasterToNumPyArray(outlet_raster, "", cols, rows, None)
nodata_o = outl_arr[0,0] # I know this is not the best way, but not an issue atm
f = numpy.argwhere(outl_arr != nodata_o)
dem_arr = arcpy.RasterToNumPyArray(inDEM, "", cols, rows, None)
# Create slope array
slope = []
x_list = numpy.arange(0,5,.01)
for i in x_list:
slope.append(-math.exp(-i))
length = len(slope)
# Make modifications
non = []
for i in range(len(f)):
index = f[i]
if dem_arr[index[0], index[1]] != nodata_dem:
r_low = index[0] - buffers
r_up = index[0] + buffers +1
c_low = index[1] - buffers
c_up = index[1] + buffers +1
data = dem_arr[r_low:r_up, c_low:c_up]
dim = math.sqrt(numpy.size(data))
n1 = math.floor(dim/2)
for j1 in range(dim):
for j2 in range(dim):
n2 = data[j1,j2]
dist = math.sqrt((n1-j1)**2 + (n1-j2)**2)*cellsz
decr = length * dist / area
if decr < length and decr != 0 and n2 != nodata_dem:
data[j1,j2] = n2 + slope[int(math.floor(decr))]
elif decr == 0:
data[j1,j2] = n2 + slope[0]
elif decr == length:
data[j1,j2] = n2 + slope[-1]
dem_arr[r_low:r_up, c_low:c_up] = data
else:
non.append(f[i])
# Trying to save the changes Alt 1.
newraster = arcpy.NumPyArrayToRaster(dem_arr, "", inDEM, "", nodata_dem)
newraster.save(outDEM)
# Trying to save Alt 2.
newdem = driver.Create(outDEM, cols, rows, 1, gdal.GDT_Float32)
newdem.SetGeoTransform(trf)
newdem.SetProjection(prj)
newband = newdem.GetRasterBand(1)
newband.WriteArray(dem_arr) #
