Python叠加两幅栅格图像

目的
现有两幅栅格图像，一个是某地区道路栅格图，一个是某地区土地利用类型图，需要将道路叠加到土地利用类型图中，即叠加后，重合的像元值以道路图为准，其余的像元值仍是土地利用类型图原有的像元值。
图1 道路信息图

图2 土地利用类型图

图3 结果图

具体实现

from gdalconst import *from osgeo import gdalimport osrimport sysimport copy#叠加两个栅格图像（一个道路栅格图，一个土地利用类型图），两幅图像重叠的像元值都是第一个图像的值，#未重叠的像元值还是土地利用类型图上的值，最终结果便是土地利用类型图上面多了道路信息。roadFile =  'E:\\Exercise\\test\\grasstest\\road_rastercalc.tif'landuseFile =  'E:\\Exercise\\test\\grasstest\\landuse.tif'roadDs = gdal.Open(roadFile, GA_ReadOnly)landuseDs = gdal.Open(landuseFile, GA_ReadOnly)if roadDs is None:    print 'Can not open ', roadFile    sys.exit(1)geotransform = roadDs.GetGeoTransform()projection=roadDs.GetProjection()cols = roadDs.RasterXSizerows = roadDs.RasterYSizeroadBand = roadDs.GetRasterBand(1)roadData = roadBand.ReadAsArray(0,0,cols,rows)roadNoData = roadBand.GetNoDataValue()landuseBand = landuseDs.GetRasterBand(1)landuseData = landuseBand.ReadAsArray(0,0,cols,rows)landuseNoData = landuseBand.GetNoDataValue()result = landuseDatafor i in range(0,rows):    for j in range(0,cols):        if(abs(roadData[i,j] - 20) < 0.0001):           result[i,j] = 20        if((abs(landuseData[i,j] - landuseNoData)>0.0001) and (abs(roadData[i,j] - roadNoData) < 0.0001)):            result[i,j] = landuseData[i,j]        if((abs(landuseData[i,j] - landuseNoData)<0.0001) and (abs(roadData[i,j] - roadNoData) < 0.0001)):            result[i,j] = landuseNoData#write result to diskresultPath = 'E:\\Exercise\\test\\grasstest\\result_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(result)    ds = Noneprint 'ok---------'