Python叠加两幅栅格图像的实现方法_Python

Python叠加两幅栅格图像的实现方法

2021-08-02 00:09hnyzwtf Python

今天小编就为大家分享一篇Python叠加两幅栅格图像的实现方法，具有很好的参考价值，希望对大家有所帮助。一起跟随小编过来看看吧

目的

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

图1 道路信息图

Python叠加两幅栅格图像的实现方法

图2 土地利用类型图

Python叠加两幅栅格图像的实现方法

图3 结果图

Python叠加两幅栅格图像的实现方法

具体实现

				?

									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---------'