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) 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, index] != nodata_dem: r_low = index - buffers r_up = index + buffers +1 c_low = index - buffers c_up = index + 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 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) #