2

My Project I am currently working at should enable cellvalue manipulation on raster files using ArcPy. My Script works fine except I am not getting the right Y-Values (northing values) for the modified raster image. I will show you my code:

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("3D")

# Local variables:
inputRaster = r'G:\Data\filegdb.gdb\ClipRaster2'
firstCoord = "33321792 6012729"
secondCoord = "33321809 6012763"
level = "5"
outputRaster = r'G:\Data\filegdb.gdb\RasterModTestData2'

# Get Raster Description:
desc = arcpy.Describe(inputRaster)
xmin = desc.extent.XMin
xmax = desc.extent.XMax
ymin = desc.extent.YMin
ymax = desc.extent.YMax
cellsize = desc.meanCellWidth
nodata = desc.noDataValue
sr = desc.SpatialReference
lowerleft = arcpy.Point(xmin,ymin) 

xwert1 = firstCoord.split(" ")
UTMRechts1 = float(xwert1[0])
UTMHoch1 = float(xwert1[1])
xwert2 = secondCoord.split(" ")
UTMRechts2 = float(xwert2[0])
UTMHoch2 = float(xwert2[1])

# Converts the Raster to a NumPyArray
raster = arcpy.Raster(inputRaster)
array = arcpy.RasterToNumPyArray(raster,nodata_to_value=raster.minimum)

# Loop through the NumPy-Array and check if the point is inside the bounding rectangle
i = ymin
countX = 0
while i < ymax:
    j = xmin
    countY = 0
    while j < xmax:

        if (j >= UTMRechts1 and j <= UTMRechts2 and i >= UTMHoch1 and i <= UTMHoch2):
            array[countX,countY] += int(level)

        countY += 1
        j += cellsize
    countX += 1
    i += cellsize

newraster = arcpy.NumPyArrayToRaster(array, lowerleft, cellsize, cellsize, nodata)
arcpy.DefineProjection_management(newraster, sr)

if arcpy.Exists(outputRaster):
    arcpy.Delete_management(outputRaster)

newraster.save(outputRaster)

As You can see, the area bounding by a rectangle (lower left: 33321792 6012729, upper right: 33321809 6012763 //UTM 33N (ETRS89) ) should be modified with raise of its elevation by 5 metres. The following screenshot shows a extract from the result raster:

Extract from result raster

Unfortunately, the lower left corner of the modified area is at (33321792 6012708 metres UTM33N/ETRS89) and the upper right corner is at (33321809 6012743 metres UTM33N/ETRS89). So the easting of the coordinate fits perfectly to the requested extent, but the northing is displaced by about 21 metres to the south. The input raster has one band, a cellsize of 0.5 m and 2254 columns and 2532 rows.

Additional Info: If I print the coordinates which are inside the bounding box, the values I want to have show up. So I guess, the mistake is in the retransformation from NumpyArray to the raster file. Has anyone experience with this problem?

Can anybody help me what to do to get the correct extent of the modified area or show me the mistake in my code?

  • did you try setting a snap raster? – mwil Oct 28 '14 at 20:44
  • you might also try setting your outputCoordinateSystem environment prior to creating it instead of defining the raster projection afterwards. – mwil Oct 28 '14 at 20:59
  • Do you mean to set the environment instead of "arcpy.DefineProjection_management(newraster, sr)" ? What do I have to set especially? – Andi90 Oct 28 '14 at 22:45
0

I solved my Problem on my own. The Problem is, that the Retransformation to the NumpyArray starts filling the array at array[0,0] on the top-left corner. So I updated my while-loop:

i = ymax
countX = 0
while i >= ymin:
    j = xmin
    countY = 0
    while j <= xmax:

    if (j >= UTMRechts1 and j <= UTMRechts2 and i >= UTMHoch1 and i <= UTMHoch2):
        array[zaehlerX,zaehlerY] += int(level)

    countY += 1
    j += cellsize
countX += 1
i -= cellsize

As you can see, I "count" from the maximum Northing-value to the minimum value. In this case, I get the right position of the modified area.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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