Python实现ARCGIS栅格计算器con函数功能

目的
现有某地区土地利用类型图，图中不同的数值代表不同的土地利用类型，如像元值为20的，我们重新赋值为0，像元值为123的，我们重新赋值为15等等。
图1 原始图像

图2 结果图像

代码如下

from gdalconst import *
from osgeo import gdal
import osr
import sys
import copy

#实现栅格计算器中的con函数功能，给图像重新赋值

landuseFile =  'E:\\Exercise\\test\\grasstest\\result_landuse.tif'
landuseDs = gdal.Open(landuseFile, GA_ReadOnly)
if landuseDs is None:
    print 'Can not open ', landuseDs
    sys.exit(1)

geotransform = landuseDs.GetGeoTransform()
projection=landuseDs.GetProjection()
cols = landuseDs.RasterXSize
rows = landuseDs.RasterYSize

landuseBand = landuseDs.GetRasterBand(1)
landuseData = landuseBand.ReadAsArray(0,0,cols,rows)
landuseNoData = landuseBand.GetNoDataValue()

resultData = landuseData

for i in range(0,rows):
    for j in range(0,cols):       
        if(abs(landuseData[i][j] - 123) < 0.0001):
            resultData[i][j] = 15
        if(abs(landuseData[i][j] - 122) < 0.0001):
           resultData[i][j] = 15
        if(abs(landuseData[i][j] - 113) < 0.0001):
           resultData[i][j] = 3
        if(abs(landuseData[i][j] - 112) < 0.0001):
           resultData[i][j] = 3
        if(abs(landuseData[i][j] - 111) < 0.0001):
           resultData[i][j] = 3
        if(abs(landuseData[i][j] - 51) < 0.0001):
           resultData[i][j] = 500
        if(abs(landuseData[i][j] - 46) < 0.0001):
           resultData[i][j] = 10
        if(abs(landuseData[i][j] - 43) < 0.0001):
           resultData[i][j] = 200
        if(abs(landuseData[i][j] - 42) < 0.0001):
           resultData[i][j] = 500
        if(abs(landuseData[i][j] - 41) < 0.0001):
           resultData[i][j] = 200
        if(abs(landuseData[i][j] - 31) < 0.0001):
           resultData[i][j] = 5
        if(abs(landuseData[i][j] - 24) < 0.0001):
           resultData[i][j] = 15
        if(abs(landuseData[i][j] - 23) < 0.0001):
           resultData[i][j] = 30
        if(abs(landuseData[i][j] - 22) < 0.0001):
           resultData[i][j] = 50
        if(abs(landuseData[i][j] - 21) < 0.0001):
           resultData[i][j] = 40
        if(abs(landuseData[i][j] - 20) < 0.0001):
           resultData[i][j] = 0

        if(abs(landuseData[i][j] - landuseNoData) < 0.0001):
            resultData[i][j] = landuseNoData

resultPath =  'E:\\Exercise\\test\\grasstest\\friction_landuse.tif'
format = "GTiff"   
driver = gdal.GetDriverByName(format)
ds = driver.Create(resultPath, cols, rows, 1, GDT_Float32)
ds.SetGeoTransform(geotransform)
ds.SetProjection(projection)
ds.GetRasterBand(1).SetNoDataValue(landuseNoData)
ds.GetRasterBand(1).WriteArray(resultData)    
ds = None
print 'ok---------'

转载自：https://blog.csdn.net/hnyzwtf/article/details/51155163