Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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) # 
share|improve this question
add comment

1 Answer

up vote 3 down vote accepted

One problem that you might be having with Alt 1 is that the second parameter in arcpy.NumpyArrayToRaster() should be an point object with the lower left coordinates of your output raster. You might try something like this:

desc = arcpy.Describe(inDEM)
xmin = desc.extent.xmin
ymin = desc.extent.ymin
cellwidth = desc.meanCellWidth
cellheight = desc.meanCellHeight
nodata = desc.noDataValue
sr = desc.SpatialReference

newraster = arcpy.NumPyArrayToRaster(dem_arr, arcpy.Point(xmin, ymin), cellwidth, cellheight, nodata)

arcpy.DefineProjection_management(newraster, sr)
newraster.save(outDEM)

Regarding Alt 2, setting your nodata value and calculating statistics might help you see your outputs.

newband.SetNoDataValue(nodata_dem)
newband.GetStatistics(0, 1)
dem = None

Hope that helps.

share|improve this answer
    
Thank you Cyrus, your suggestions on Alt 1 worked perfectly! Now the new raster is located where it should be. –  Martin Nov 8 '12 at 8:21
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.