目的
现有某地区土地利用类型图,图中不同的数值代表不同的土地利用类型,如像元值为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---------'