Python叠加两幅栅格图像的实现方法
作者:hnyzwtf 时间:2023-07-25 17:14:06
目的
现有两幅栅格图像,一个是某地区道路栅格图,一个是某地区土地利用类型图,需要将道路叠加到土地利用类型图中,即叠加后,重合的像元值以道路图为准,其余的像元值仍是土地利用类型图原有的像元值。
图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---------'
来源:https://blog.csdn.net/hnyzwtf/article/details/51155090
标签:Python,叠加,栅格,图像
![](/images/zang.png)
![](/images/jiucuo.png)
猜你喜欢
Python实现的列表排序、反转操作示例
2023-06-19 11:14:27
Python统计词频并绘制图片(附完整代码)
2022-01-01 09:28:24
![](https://img.aspxhome.com/file/2023/7/110447_0s.jpg)
Go调度器学习之goroutine调度详解
2024-04-30 10:06:10
![](https://img.aspxhome.com/file/2023/4/130694_0s.jpg)
Python基于locals返回作用域字典
2021-05-17 22:02:43
Python实现爬取网页中动态加载的数据
2021-08-11 18:35:27
![](https://img.aspxhome.com/file/2023/7/102857_0s.jpg)
asp中数组的用法
2008-05-12 22:29:00
浅谈keras中的Merge层(实现层的相加、相减、相乘实例)
2021-04-03 15:30:48
关注oblog 关注xml-rpc 自己开发客户端小记
2009-10-04 20:27:00
thinkphp5框架调用其它控制器方法 实现自定义跳转界面功能示例
2024-05-03 15:52:31
使用Python的Django和layim实现即时通讯的方法
2022-05-29 00:41:33
![](https://img.aspxhome.com/file/2023/0/77570_0s.jpg)
AI与Python人工智能启发式搜索概念理解
2021-08-05 03:56:28
MySQL 8 新特性之Invisible Indexes
2024-01-19 04:56:06
浅谈Python批处理文件夹中的txt文件
2023-03-17 22:58:28
![](https://img.aspxhome.com/file/2023/3/65033_0s.png)
详解MYSQL的备份还原(PHP实现)
2023-10-06 10:44:23
Python创建xml的方法
2021-11-25 05:59:44
Python使用pyecharts绘制世界地图,省级地图,城市地图实例详解
2022-08-08 06:37:22
![](https://img.aspxhome.com/file/2023/9/72119_0s.jpg)
python excel多行合并的方法
2021-09-25 22:36:29
![](https://img.aspxhome.com/file/2023/3/81823_0s.png)
webpack cjs运行时分析示例详解
2024-04-19 09:51:56
![](https://img.aspxhome.com/file/2023/2/136072_0s.jpg)
Go 语言的 :=的具体使用
2024-05-10 13:58:56
![](https://img.aspxhome.com/file/2023/9/127329_0s.png)
Python可视化程序调用流程解析
2022-07-18 15:53:07
![](https://img.aspxhome.com/file/2023/5/133235_0s.jpg)