Python 实现遥感影像波段组合的示例代码

最近要做个遥感相关的小系统，需要波段组合功能，网上找了可以使用ArcGIS安装时自带的arcpy包，但是Python3.7不能使用现有ArcGIS10.2版本，也不想再装其他版本，所以只能自己想了个办法解决。不过有点笨啊。

思路是：

1.读取需要组合遥感影像波段（此处用OLI）

2.创建数组，把读取的波段按序放进去

3.写入文件，写成tif多波段数据

上代码：

from osgeo import gdal
import os
import numpy as np

class GRID:

  #读图像文件
  def read_img(self,filename):
    dataset=gdal.Open(filename)    #打开文件

    im_width = dataset.RasterXSize  #栅格矩阵的列数
    im_height = dataset.RasterYSize  #栅格矩阵的行数

    im_geotrans = dataset.GetGeoTransform() #仿射矩阵
    im_proj = dataset.GetProjection() #地图投影信息
    im_data = dataset.ReadAsArray(0,0,im_width,im_height) #将数据写成数组，对应栅格矩阵

    del dataset #关闭对象，文件dataset
    return im_proj,im_geotrans,im_data,im_width,im_height

  #写文件，以写成tif为例
  def write_img(self,filename,im_proj,im_geotrans,im_data):

    #判断栅格数据的数据类型
    if 'int8' in im_data.dtype.name:
      datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
      datatype = gdal.GDT_UInt16
    else:
      datatype = gdal.GDT_Float32

    #判读数组维数
    if len(im_data.shape) == 3:
      im_bands, im_height, im_width = im_data.shape
    else:
      im_bands, (im_height, im_width) = 1,im_data.shape

    #创建文件
    driver = gdal.GetDriverByName("GTiff")  #数据类型必须有，因为要计算需要多大内存空间
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)       #写入仿射变换参数
    dataset.SetProjection(im_proj)          #写入投影

    if im_bands == 1:
      dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据
    else:
      for i in range(im_bands):
        dataset.GetRasterBand(i+1).WriteArray(im_data[i])

    del dataset

if __name__ == "__main__":
  os.chdir(r'E:\Python\temp\data')            #切换路径到待处理图像所在文件夹
  run = GRID()
  #第一步
  proj,geotrans,data1,row1,column1 = run.read_img('Band_5_Clip.tif') #读数据
  proj,geotrans,data2,row2,column2= run.read_img('Band_4_Clip.tif') # 读数据
  proj,geotrans,data3,row3,column3 = run.read_img('Band_3_Clip.tif') # 读数据

  #第二步
  data = np.array((data1, data2, data3),dtype = data1.dtype) #按序将3个波段像元值放入

  #第三步
  run.write_img('com543.tif', proj, geotrans, data) # 写为3波段数据

OK！！和ArcGIS处理的对比一下，发现差别不大(上：ArcMap   下：Python）。

方法较笨，如果各位大神有更好的方法，我们可以私下交流交流。

以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持我们。

赞 (0)

相关推荐

  • 利用python GDAL库读写geotiff格式的遥感影像方法

    如下所示: from osgeo import gdal import numpy as np def read_tiff(inpath): ds=gdal.Open(inpath) row=ds.RasterXSize col=ds.RasterYSize band=ds.RasterCount geoTransform=ds.GetTransform() proj=ds.GetTransform() data=np.zeros([row,col,band]) for i in range(b

  • Python中使用OpenCV库来进行简单的气象学遥感影像计算

    Python中使用OpenCV库来进行简单的气象学遥感影像计算

    OpenCV的全称是Open Source Computer Vision Library,是一个跨平台的计算机视觉库.OpenCV是由英特尔公司发起并参与开发,以BSD许可证授权发行,可以在商业和研究领域中免费使用.OpenCV可用于开发实时的图像处理.计算机视觉以及模式识别程序.该程序库也可以使用英特尔公司的IPP进行加速处理. OpenCV用C++语言编写,它的主要接口也是C++语言,但是依然保留了大量的C语言接口.该库也有大量的Python, Java and MATLAB/OCTAVE

  • Python 实现遥感影像波段组合的示例代码

    最近要做个遥感相关的小系统,需要波段组合功能,网上找了可以使用ArcGIS安装时自带的arcpy包,但是Python3.7不能使用现有ArcGIS10.2版本,也不想再装其他版本,所以只能自己想了个办法解决.不过有点笨啊. 思路是: 1.读取需要组合遥感影像波段(此处用OLI) 2.创建数组,把读取的波段按序放进去 3.写入文件,写成tif多波段数据 上代码: from osgeo import gdal import os import numpy as np class GRID: #读图像

  • 详解Python修复遥感影像条带的两种方式

    GDAL修复Landsat ETM+影像条带 Landsat7 ETM+卫星影像由于卫星传感器故障,导致此后获取的影像出现了条带.如下图所示, 影像中均匀的布满条带. 使用GDAL修复影像条带的代码如下: def gdal_repair(tif_name, out_name, bands): """ tif_name(string): 源影像名 out_name(string): 输出影像名 bands(integer): 影像波段数 """ #

  • Python实现清理重复文件功能的示例代码

    Python实现清理重复文件功能的示例代码

    目录 前置 查找.删除重复文件 GUI制作 GUI界面设计 逻辑设计 效果展示 在电脑上或多或少的存在一些重复文件,体积小的倒没什么,如果体积大的就很占内存了,而如果自己一个一个查看文件是否重复,然后再删除,还是很要命的. 为此,我用python制作了一个删除重复文件的小工具,核心代码很简单,就十行代码,不管什么类型的文件都可以一键删除! 前置 PySimpleGUI库用来创建可视化界面,os操作文件,只需要这两个库: import os import PySimpleGUI as sg os为

  • Python实现邮件的批量发送的示例代码

    1 发送文本信息 '''加密发送文本邮件''' def sendEmail(from_addr,password,to_addr,smtp_server): try: msg = MIMEText('你好,来自信息化工程所的问候...', 'plain', 'utf-8') # 文本邮件 # msg = MIMEText('<html><body><h1>你好</h1>' + '<p>send by <a href="http:/

  • python画双y轴图像的示例代码

    很多时候可能需要在一个图中画出多条函数图像,但是可能y轴的物理含义不一样,或是数值范围相差较大,此时就需要双y轴. matplotlib和seaborn都可以画双y轴图像. 一个例子: import seaborn as sns import matplotlib.pyplot as plt # ax1 for KDE, ax2 for CDF f, ax1 = plt.subplots() ax1.grid(True) # ax1.set_ylim(0, 1) ax1.set_ylabel('

  • 用python实现刷点击率的示例代码

    背景 同事的老爸参加微信的一个活动,需要刷点击率,因此,写了一个程序助之. 准备 微信活动也是有真实地址的. 通过mitmproxy(man in the middle proxy)的方式,可以获取微信获取网页的真实地址(url). 完整可运行代码 import os import time import argparse import platform def visit_win(url, times, duration): import urllib2 def _visit_win(): t

  • Python Learning 列表的更多操作及示例代码

    遍历列表-for循环 列表中存储的元素可能非常多,如果想一个一个的访问列表中的元素,可能是一件十分头疼的事.那有没有什么好的办法呢?当然有!使用 for循环 假如有一个食物名单列表,通过 for循环 将列表中的食物名称都打印出来 # 定义一个食物名单列表 foods = ['potato', 'tomato', 'noodles', 'apple', 'pizza'] # 循环访问foods列表 for food in foods: print(food) 输出: potato  tomato 

  • Python实现RabbitMQ6种消息模型的示例代码

    RabbitMQ与Redis对比 ​ RabbitMQ是一种比较流行的消息中间件,之前我一直使用redis作为消息中间件,但是生产环境比较推荐RabbitMQ来替代Redis,所以我去查询了一些RabbitMQ的资料.相比于Redis,RabbitMQ优点很多,比如: 具有消息消费确认机制 队列,消息,都可以选择是否持久化,粒度更小.更灵活. 可以实现负载均衡 RabbitMQ应用场景 异步处理:比如用户注册时的确认邮件.短信等交由rabbitMQ进行异步处理 应用解耦:比如收发消息双方可以使用

  • Python爬虫实现vip电影下载的示例代码

    爬虫目的 实现对各大视频网站vip电影的下载,因为第三方解析网站并没有提供下载的渠道,因此想要实现电影的下载. 实现思路 1.选择一个合适的vip解析网站,这里选择了无名小站的接口,因为尝试了很多网站,有些网站想要爬取很困难,无名小站相对简单,接口为www.wmxz.wang/video.php?url=[vip电影的链接] 2.利用Fiddler进行抓包,模拟浏览器发送post请求,获取电影实际下载地址. 3.使用PyQt5进行包装,实现多样化的功能.(可选) 页面分析 我使用Fiddler抓

  • Python实现自动打开电脑应用的示例代码

    由于时间原因,有时候可能会错过某个上网课的时间段.因此想要实现自动定时启动DingDing. 新手一枚,如有不当勿喷望大佬指正. 自动打开DingDing可以由两种方法实现: 通过找出找出软件在电脑中快捷方式的位置(电脑屏幕中的坐标),使用代码模拟鼠标进行双击打开. 通过输入软件在电脑中的安装路径打开软件. 1.第一种方法: ​在python中,使用pip install pyautogui 安装第三方库,在此库中,可以使用pyautogui.position()获取当前鼠标放置位置的坐标.我们

随机推荐