现有两幅栅格图像,一个是某地区道路栅格图,一个是某地区土地利用类型图,需要将道路叠加到土地利用类型图中,即叠加后,重合的像元值以道路图为准,其余的像元值仍是土地利用类型图原有的像元值。
<强>图1道路信息图
<强>图2土地利用类型图
<强>图3结果图
<强>具体实现
从gdalconst进口*
从osgeo进口gdal
进口osr
导入系统
进口复制
#叠加两个栅格图像(一个道路栅格图,一个土地利用类型图),两幅图像重叠的像元值都是第一个图像的值,
#未重叠的像元值还是土地利用类型图上的值,最终结果便是土地利用类型图上面多了道路信息。
roadFile=' E: \ \ \ \运动测试\ \ grasstest \ \ road_rastercalc.tif '
landuseFile=' E: \ \ \ \运动测试\ \ grasstest \ \ landuse.tif '
roadDs=gdal。Open (roadFile GA_ReadOnly)
landuseDs=gdal。Open (landuseFile GA_ReadOnly)
如果roadDs没有:
打印“无法打开”,roadFile
sys.exit (1)
geotransform=roadDs.GetGeoTransform ()
投影=roadDs.GetProjection ()
关口=roadDs.RasterXSize
行=roadDs.RasterYSize
roadBand=roadDs.GetRasterBand (1)
roadData=https://www.yisu.com/zixun/roadBand.ReadAsArray(0, 0,关口,行)
roadNoData=roadBand.GetNoDataValue ()
landuseBand=landuseDs.GetRasterBand (1)
landuseData=landuseBand.ReadAsArray(0, 0,关口,行)
landuseNoData=landuseBand.GetNoDataValue ()
结果=landuseData
因为我在范围(0,行):
j的范围(0,峡路):
如果(abs (roadData (i, j) - 20) <0.0001):
结果(i, j)=20
如果((abs (landuseData (i, j) - landuseNoData)> 0.0001)和(abs (roadData (i, j) - roadNoData) <0.0001)):
结果(i, j)=landuseData (i, j)
如果((abs (landuseData (i, j) - landuseNoData) <0.0001)和(abs (roadData (i, j) - roadNoData) <0.0001)):
结果(i, j)=landuseNoData
#结果写入到磁盘
resultPath=' E: \ \ \ \运动测试\ \ grasstest \ \ result_landuse.tif '
格式=" GTiff "
司机=gdal.GetDriverByName(格式)
ds=司机。创建(resultPath关口,行,1,GDT_Float32)
ds.SetGeoTransform (geotransform)
ds.SetProjection(投影)
ds.GetRasterBand (1) .SetNoDataValue (landuseNoData)
ds.GetRasterBand (1) .WriteArray(结果)
ds=没有
打印“好吧- - - - - - - - - - - -”
之前
以上这篇Python叠加两幅栅格图像的实现方法就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持。
