python图像加网格_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叠加两幅栅格图像的实现方法

本文地址: http://www.cppcns.com/jiaoben/python/264714.html