Landsat8处理小工具(python)

再更新

决定将重构的代码更新在这篇博客中。现在实现的功能：读、写、多波段的辐射定标（可定表观辐射亮度或表观反射率）、图像分片功能（全部带坐标）。
读取可选多光谱波段、全色波段、热红外波段，读入numpy数组。
辐射定标后读入numpy数组。

更新

决定不再更新这段代码，而是重构代码。因为强行使用装饰器等，工具使用起来十分蹩脚。至于新的重构代码，功能会更加齐全，会更新在新的博客中。

之前内容

暂时有了读、写、分块（写出时暂无坐标）、多光谱辐射定标等功能。
下面会不停更新此博客，进而更新代码。

至于为什么用装饰器等，主要是为了练习使用而已。

代码

# -*-coding:utf-8-*-
#
# @author:rsstudent
#import numpy as np
from osgeo import gdal
from osgeo import gdal_array
import cv2 as cv
import math
from numba import jit
from functools import wraps
import osclass Landsat8Tools(object):def __init__(self):self.LANDSAT8_MULTI_BAND_NUM = 7self.LANDSAT8_THERMAL_BAND_NUM = 2self.LANDSAT8_PAN_BAND_NUM = 1self.LANDSAT8_PAN_BAND = 8self.LANDSAT8_THERMAL_START_BAND = 10self.multi_nan_position = []self.pan_nan_position = []self.thermal_nan_position = []self.process_image_nan_position = []def __get_base_name(self, full_name):base_names = full_name.split("_")[:-1]base_name = ""for string in base_names:base_name += (string + "_")return base_namedef read_multi_band_to_image(self, mtl_file_path):multi_band_file_names = []base_name = self.__get_base_name(mtl_file_path)for band in range(self.LANDSAT8_MULTI_BAND_NUM):band_name = base_name + "B" + str(band + 1) + ".tif"multi_band_file_names.append(band_name)if len(multi_band_file_names) == 0:raise Exception("The filename is incorrect!")dataset = gdal.Open(multi_band_file_names[0], gdal.GA_ReadOnly)if dataset == None:raise Exception("Fail to open the image file.")height = dataset.RasterYSizewidth = dataset.RasterXSizeimage = np.zeros((height, width, self.LANDSAT8_MULTI_BAND_NUM), dtype = np.float16)projection = dataset.GetProjection()geotransform = dataset.GetGeoTransform()del datasetprint("read multi band start, please wait...")for band in range(self.LANDSAT8_MULTI_BAND_NUM):ds = gdal.Open(multi_band_file_names[band], gdal.GA_ReadOnly)band_image = ds.GetRasterBand(1)image[:, :, band] = band_image.ReadAsArray()del dsprint("read  multi band finish.")self.multi_nan_position = np.where(image == 0)image[self.multi_nan_position] = np.nanreturn image, projection, geotransformdef read_pan_band_to_image(self, mtl_file_path):base_name = self.__get_base_name(mtl_file_path)pan_file_name = base_name + "B8" + ".tif"dataset = gdal.Open(pan_file_name, gdal.GA_ReadOnly)if dataset == None:raise Exception("Image name error.")else:datatype = np.float16height = dataset.RasterYSizewidth = dataset.RasterXSizeprojection = dataset.GetProjection()geotransform = dataset.GetGeoTransform()pan_image = np.zeros((height, width, 1), dtype = datatype)print("read pan start, please wait...")band_data = dataset.GetRasterBand(1)pan_image[:, :, 0] = band_data.ReadAsArray()del datasetself.pan_nan_position = np.where(pan_image == 0)pan_image[self.pan_nan_position] = np.nanprint("read pan finish.")return pan_image, projection, geotransformdef read_thermal_band_to_image(self, mtl_file_path):nan_position = []base_name = self.__get_base_name(mtl_file_path)thermal_band10_file_name = base_name + "B10" + ".tif"thermal_band11_file_name = base_name + "B11" + ".tif"names = []names.append(thermal_band10_file_name)names.append(thermal_band11_file_name)dataset = gdal.Open(thermal_band10_file_name, gdal.GA_ReadOnly)if dataset == None:raise Exception("Image name error.")else:datatype = np.float16height = dataset.RasterYSizewidth = dataset.RasterXSizeprojection = dataset.GetProjection()geotransform = dataset.GetGeoTransform()thermal_image = np.zeros((height, width, 2), dtype = datatype)del datasetprint("read thermal band start, please wait...")for band in range(self.LANDSAT8_THERMAL_BAND_NUM):ds = gdal.Open(names[band], gdal.GA_ReadOnly)band_data = ds.GetRasterBand(1)thermal_image[:, :, band] = band_data.ReadAsArray()del dsself.thermal_nan_position = np.where(thermal_image == 0)thermal_image[self.thermal_nan_position] == np.nanprint("read thermal band finish.", end = "")return thermal_image, projection, geotransformdef read_single_band_to_image(self, mtl_file_path, band_num):base_name =self.__get_base_name(mtl_file_path)nan_position = []fullbandname = base_name + "B" + str(band_num) + ".tif"dataset = gdal.Open(fullbandname, gdal.GA_ReadOnly)if dataset == None:raise Exception("Image name error.")else:datatype = np.float16height = dataset.RasterYSizewidth = dataset.RasterXSizeprojection = dataset.GetProjection()geotransform = dataset.GetGeoTransform()band_image = np.zeros((height, width, 1), dtype = datatype)print("read single band start")band_data = dataset.GetRasterBand(1)band_image[:, :, 0] = band_data.ReadAsArray()del datasetprint("read single band finish.")self.process_image_nan_position = np.where(band_image == 0)band_image[self.process_image_nan_position] == np.nanreturn band_image, projection, geotransformdef save(self, save_path, image, projection, geotransform, format = 'GTiff'):datatype = gdal.GDT_Float32DIMENSION_OF_IMAGE = 3if len(image.shape) != DIMENSION_OF_IMAGE:raise Exception("The dimension of the image is incorrect.")else:height = image.shape[0]width = image.shape[1]channels = image.shape[2]driver = gdal.GetDriverByName(format)ds_to_save = driver.Create(save_path, width, height, channels, datatype)ds_to_save.SetGeoTransform(geotransform)ds_to_save.SetProjection(projection)print("save tool start, please wait...")for band in range(channels):ds_to_save.GetRasterBand(band + 1).WriteArray(image[:, :, band])ds_to_save.FlushCache()print("save finish.")del imagedel ds_to_savedef radiometric_calibration(self, mtl_file_path, save_folder, cali_type = "radiance"):f = open(mtl_file_path, 'r')metadata = f.readlines()f.close()radiance_multi_paras = []radiance_add_paras = []reflectance_multi_paras = []reflectance_add_paras = []radiance_paras_start_line = 0reflectance_paras_start_line = 0for lines in metadata:test_line = lines.split("=")if test_line[0] == '    REFLECTANCE_MULT_BAND_1 ':breakelse:reflectance_paras_start_line += 1for lines in range(reflectance_paras_start_line, reflectance_paras_start_line + 9):parameter = float(metadata[lines].split('=')[1])reflectance_multi_paras.append(parameter)for lines in range(reflectance_paras_start_line + 9, reflectance_paras_start_line + 18):parameter = float(metadata[lines].split('=')[1])reflectance_add_paras.append(parameter)        for lines in metadata:test_line = lines.split("=")if test_line[0] == '    RADIANCE_MULT_BAND_1 ':breakelse:radiance_paras_start_line += 1for lines in range(radiance_paras_start_line, radiance_paras_start_line + 11):parameter = float(metadata[lines].split('=')[1])radiance_multi_paras.append(parameter)for lines in range(radiance_paras_start_line + 11, radiance_paras_start_line + 22):parameter = float(metadata[lines].split('=')[1])radiance_add_paras.append(parameter)if cali_type != "radiance" and cali_type != "reflectance":raise Exception("cali_type is incorrect.")if cali_type == "radiance":multi_image, multi_projection, multi_geotransform = self.read_multi_band_to_image(mtl_file_path)print("radiance radiomereic calibration start.")for band in range(self.LANDSAT8_MULTI_BAND_NUM):gain = radiance_multi_paras[band]offset = radiance_multi_paras[band]multi_image[:, :, band] = multi_image[:, :, band]*gain + offsetmulti_image[self.multi_nan_position] = np.nanmulti_image_name = "radiance_multi.tif"save_path = os.path.join(save_folder, multi_image_name)self.save(save_path, multi_image, multi_projection, multi_geotransform)del multi_imagepan_image, pan_projection, pan_geotransform = self.read_pan_band_to_image(mtl_file_path)pan_gain = radiance_multi_paras[7]pan_offset = radiance_add_paras[7]pan_image = pan_image * gain + offsetpan_image[self.pan_nan_position] = np.nanpan_image_name = "radiance_pan.tif"save_path = os.path.join(save_folder, pan_image_name)self.save(save_path, pan_image, pan_projection, pan_geotransform)del pan_imagethermal_image, thermal_projection, thermal_geotransform = self.read_thermal_band_to_image(mtl_file_path)thermal_multi_paras = []thermal_add_paras = []thermal_multi_paras.append(radiance_multi_paras[9])thermal_multi_paras.append(radiance_multi_paras[10])thermal_add_paras.append(radiance_add_paras[9])thermal_add_paras.append(radiance_add_paras[10])for band in range(self.LANDSAT8_THERMAL_BAND_NUM):thermal_image[:, :, band] = thermal_image[:, :, band] * thermal_multi_paras[band] + thermal_add_paras[band]thermal_image[self.thermal_nan_position] = np.nanthermal_image_name = "radiance_thermal.tif"save_path = os.path.join(save_folder, thermal_image_name)self.save(save_path, thermal_image, thermal_projection, thermal_geotransform)print("ridiometric calibraion finish.")del thermal_imageelse:multi_image, multi_projection, multi_geotransform = self.read_multi_band_to_image(mtl_file_path)print("reflectance radiomereic calibration start.")for band in range(self.LANDSAT8_MULTI_BAND_NUM):gain = reflectance_multi_paras[band]offset = reflectance_multi_paras[band]multi_image[:, :, band] = multi_image[:, :, band]*gain + offsetmulti_image[self.multi_nan_position] = np.nanmulti_image_name = "reflectance_multi.tif"save_path = os.path.join(save_folder, multi_image_name)self.save(save_path, multi_image, multi_projection, multi_geotransform)del multi_imagepan_image, pan_projection, pan_geotransform = self.read_pan_band_to_image(mtl_file_path)pan_gain = reflectance_multi_paras[7]pan_offset = reflectance_add_paras[7]pan_image = pan_image * gain + offsetpan_image[self.pan_nan_position] = np.nanpan_image_name = "reflectance_pan.tif"save_path = os.path.join(save_folder, pan_image_name)self.save(save_path, pan_image, pan_projection, pan_geotransform)del pan_imagethermal_image, thermal_projection, thermal_geotransform = self.read_thermal_band_to_image(mtl_file_path)thermal_multi_paras = []thermal_add_paras = []thermal_multi_paras.append(radiance_multi_paras[9])thermal_multi_paras.append(radiance_multi_paras[10])thermal_add_paras.append(radiance_add_paras[9])thermal_add_paras.append(radiance_add_paras[10])for band in range(self.LANDSAT8_THERMAL_BAND_NUM):thermal_image[:, :, band] = thermal_image[:, :, band] * thermal_multi_paras[band] + thermal_add_paras[band]thermal_image[self.thermal_nan_position] = np.nanthermal_image_name = "radiance_thermal.tif"save_path = os.path.join(save_folder, thermal_image_name)self.save(save_path, thermal_image, thermal_projection, thermal_geotransform)print("radiometric calibraion finish.")del thermal_imagedef cut_to_tile_with_geoinfo(self, save_folder, image_size, image, projection, geotransform, format = "GTiff"):print("cut to tile with geoinfo tool start.")DIMENSION_OF_IMAGE = 3if len(image_size) !=  DIMENSION_OF_IMAGE:raise Exception("your image dimension is incorrect.Need 3-d array.")image_height = image.shape[0]image_width = image.shape[1]image_channels = image.shape[2]tile_height = image_size[0]tile_width = image_size[1]channels = image_size[2]if image_channels != channels:raise Exception("Your image channels isn't match the required channels.")num_rows = math.floor(image_height/tile_height)num_cols = math.floor(image_width/tile_width)useful_height = num_rows * tile_heightuseful_width = num_cols * tile_widthimage = np.array(image[:useful_height, :useful_width, :])x_origin_point = geotransform[0]x_pixel_size = geotransform[1]x_ro = geotransform[2]y_origin_point = geotransform[3]y_ro = geotransform[4]y_pixel_size = geotransform[5]i = 0x_offset = 0y_offset = 0for tile_row in range(num_rows):for tile_col in range(num_cols):tile = image[tile_row*tile_width: (tile_row+1)*tile_width, tile_col*tile_height:(tile_col+1)*tile_height, :]filename = str(i) + '.tif'path = os.path.join(save_folder, filename)tile = np.float32(tile)tile_x_origin_point = x_origin_point + x_pixel_size * tile_col * tile_widthtile_y_origin_point = y_origin_point + y_pixel_size * tile_row * tile_heighttile_geotransform = (tile_x_origin_point, x_pixel_size, x_ro, tile_y_origin_point, y_ro, y_pixel_size)self.save(path, tile, projection, tile_geotransform, format)i += 1del tileprint("cut to tile with geoinfo tool finish.")del imageif __name__ == "__main__":tool = Landsat8Tools()mtl_file_path = "F:\SRGAN_program\dataset\LC81290352019095LGN00\LC08_L1TP_129035_20190405_20190422_01_T1_MTL.txt"save_folder = "F:\\SRGAN_program\\dataset\\LC81290352019095LGN00\\tiles"save_folder1 = "F:\\SRGAN_program\\dataset\\LC81290352019095LGN00\\test\\ppan"pan, projection, geotransform = tool.read_pan_band_to_image(mtl_file_path)tool.cut_to_tile_with_geoinfo(save_folder, (256, 256, 1), pan, projection, geotransform)#tool.save(save_folder1, pan, projection, geotransform)

Landsat8处理小工具(python)相关推荐

pdf合并小工具python
这里写自定义目录标题前言代码前言使用python写了个合并pdf的小工具,共享代码相对比较简陋,将就着看吧,需要的可以自己再润色一下代码 from os import rename fro ...
ADB操作手机的一个界面小工具(python实现)
我们经常使用adb命令操作手机,因此我突然想到做一个界面把这些命令用界面点击的形式操作,这样可以简化我们平时敲命令的时间,工具的功能尚不完善,先总结一下思路.首先先把工具的界面展示一下: 首先讲一下工 ...
python编写测试小工具-Python与游戏测试（小工具篇）
import datetime import paramiko import time import os class ScanError(object): def __init__(self): s ...
炉石战棋一键断网小工具(python)
通过pyqt5编写的界面,通过os对防火墙入站规则进行开启与关闭注意事项: 1.对程序路径进行符号转换,否则失效 2.关闭相关杀毒软件测试 3.以管理员身份运行,打包文件exe时同理开启管理员模式 ...
python 写入excel_实用小工具python数组快速写入excel表格
在进行数值计算时有时会产生大量的数据,有时候需要与其他软件计算的数据进行对比,这时候将数据输出excel表格就很重要了. 借助openpyxl可以方便地输出相应的数据,精简版如下: from open ...
python自动翻译小工具_Python实现翻译小工具
一.背景利用Requests模块获取有道词典web页面的post信息,BeautifulSoup来获取需要的内容,通过tkinter模块生成gui界面. 二.代码 git源码地址 Python实现翻 ...
把python语言翻译出来_Python语言实现翻译小工具（Python打包成exe文件）
本文主要向大家介绍了Python语言实现翻译小工具(Python打包成exe文件),通过具体的内容向大家展示,希望对大家学习Python语言有所帮助. 1.环境 windows10 python3.5 ...
python提高办公效率-几个可以提高工作效率的Python内置小工具
在这篇文章里,我们将会介绍4个Python解释器自身提供的小工具.这些小工具在笔者的日常工作中经常用到,减少了各种时间的浪费,然而,却很容易被大家忽略.每当有新来的同事看到我这么使用时,都忍不住感叹, ...
python对工作效率的提升_使用了这个几个Python内置小工具，可以让你的工作效率提升一倍...
使用了这个几个Python内置小工具,可以让你的工作效率提升一倍我们将会详情4个Python解释器自身提供的小工具. 这些小工具在笔者的日常工作中经常使用到, 减少了各种时间的白费, 然而,却很容易 ...
自己整理实现的python小工具
文章目录记录一些自己整理实现的python小工具 python获取文件路径 pytho使用opencv进行图像拼接记录一些自己整理实现的python小工具 python获取文件路径因为有的程序需 ...

Landsat8处理小工具(python)

再更新

更新

之前内容

代码

Landsat8处理小工具(python)相关推荐

最新文章

热门文章