@本文来源于公众号:csdn2299,喜欢可以关注公众号 程序员学府
今天小编就为大家分享一篇Python叠加两幅栅格图像的实现方法,具有很好的参考价值,希望对大家有所帮助。一起跟随小编过来看看吧
目的
现有两幅栅格图像,一个是某地区道路栅格图,一个是某地区土地利用类型图,需要将道路叠加到土地利用类型图中,即叠加后,重合的像元值以道路图为准,其余的像元值仍是土地利用类型图原有的像元值。
图1 道路信息图
图2 土地利用类型图
图3 结果图
具体实现
from gdalconst import *
from osgeo import gdal
import osr
import sys
import 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.RasterXSize
rows = roadDs.RasterYSize
roadBand = 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 = landuseData
for 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 disk
resultPath = '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 = None
print 'ok---------'
非常感谢你的阅读
大学的时候选择了自学python,工作了发现吃了计算机基础不好的亏,学历不行这是没办法的事,只能后天弥补,于是在编码之外开启了自己的逆袭之路,不断的学习python核心知识,深入的研习计算机基础知识,整理好了,我放在我们的微信公众号《程序员学府》,如果你也不甘平庸,那就与我一起在编码之外,不断成长吧!
其实这里不仅有技术,更有那些技术之外的东西,比如,如何做一个精致的程序员,而不是“屌丝”,程序员本身就是高贵的一种存在啊,难道不是吗?[点击加入]
想做你自己想成为高尚人,加油!