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) UTMHoch1 = float(xwert1) xwert2 = secondCoord.split(" ") UTMRechts2 = float(xwert2) UTMHoch2 = float(xwert2) # 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:
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?