程序员经验分享-hwhale

Python遥感开发之分段读取和保存遥感数据

2023-11-15


前言：当遇到批量读取大量遥感数据进行运算的时候，如果不进行分段读取操作的话，电脑内存可能面临着不够使用的情况，所以我们要进行分段读取数据然后进行运算，运算结束之后把这段数据保存成tif文件。

1 分段读取数据

在这里插入图片描述
如图所示，有三个这样的数据，且该数据为5600行6800列，我们可以分成10个批次分段读取该TIF数据，10个批次以此为0,560,1120,1680,2240,2800,3360,3920,4480,5040,5600。
代码实现：

import os
import numpy as np
from osgeo import gdal, gdalnumeric

def read_tif(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize#图像长度
    row = dataset.RasterYSize#图像宽度
    geotrans = dataset.GetGeoTransform()#读取仿射变换
    proj = dataset.GetProjection()#读取投影
    data = dataset.ReadAsArray()#转为numpy格式
    data = data.astype(np.float32)
    # a = data[0][0]
    # data[data == a] = np.nan
    data[data <= -3] = np.nan
    return [col, row, geotrans, proj, data]

def read_tif02(file):
    data = gdalnumeric.LoadFile(file)
    data = data.astype(np.float32)
    # a = data[0][0]
    # data[data == a] = np.nan
    data[data <= -3] = np.nan
    return data

def get_all_file_name(ndvi_file):
    list1=[]
    if os.path.isdir(ndvi_file):
        fileList = os.listdir(ndvi_file)
        for f in fileList:
            file_name= ndvi_file+"\\"+f
            list1.append(file_name)
        return list1
    else:
        return []

if __name__ == '__main__':
    file_ndvi = r"D:\AAWORK\work\研究方向\研究方向01作物分类内容\风云数据\MERSI-II植被指数旬产品(250M)\NDVI主合成"
    file_out = r"D:\AAWORK\work\2021NDVISUM.tif"
    ndvi_file_list = get_all_file_name(file_ndvi)
    col00, row00, geotrans00, proj00, data_00_ndvi = read_tif(ndvi_file_list[0])
    data_01_ndvi = read_tif02(ndvi_file_list[1])
    data_02_ndvi = read_tif02(ndvi_file_list[2])
    list_row = [0,560,1120,1680,2240,2800,3360,3920,4480,5040,5600]
    for index,i in enumerate(list_row):
        if index <= len(list_row)-2:
            print(list_row[index],list_row[index+1])
            #分段进行操作
            # sum_list = get_section(list_row[index],list_row[index+1],col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)
            # 分段进行保存
            # save_section(sum_list, list_row[index], list_row[index+1], col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)

在这里插入图片描述

2 实现分批读取数据以及进行计算

拿到开始的行和结束的行数，进行分批读取数据并进行计算，（这里求和求的是整数，如有需要的话可以自己更改）代码如下：

import os
import tensorflow as tf
import numpy as np
import pandas as pd
from osgeo import gdal, gdalnumeric

def get_sum_list(data_list):
    list1 = []
    for data in data_list:
        sum = 0
        for d in data:
            if not np.isnan(d):
                sum = sum+d
        list1.append(int(sum))
    return list1

def read_tif(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize#图像长度
    row = dataset.RasterYSize#图像宽度
    geotrans = dataset.GetGeoTransform()#读取仿射变换
    proj = dataset.GetProjection()#读取投影
    data = dataset.ReadAsArray()#转为numpy格式
    data = data.astype(np.float32)
    # a = data[0][0]
    # data[data == a] = np.nan
    data[data <= -3] = np.nan
    return [col, row, geotrans, proj, data]

def read_tif02(file):
    data = gdalnumeric.LoadFile(file)
    data = data.astype(np.float32)
    # a = data[0][0]
    # data[data == a] = np.nan
    data[data <= -3] = np.nan
    return data

def get_all_file_name(ndvi_file):
    list1=[]
    if os.path.isdir(ndvi_file):
        fileList = os.listdir(ndvi_file)
        for f in fileList:
            file_name= ndvi_file+"\\"+f
            list1.append(file_name)
        return list1
    else:
        return []

def get_nan_sum(ndvi_list):
    """
    得到NAN数据的个数
    :param ndvi_list:
    :return:
    """
    count = 0
    for ndvi in ndvi_list:
        if np.isnan(ndvi):
            count = count+1
    return count

def get_section(row0, row1, col1,data1,data2,data3):
    """
    分段读取数据，读取的数据进行计算
    :param row0:
    :param row1:
    :param col1:
    :param data1:
    :param data2:
    :param data3:
    :return:
    """
    list1 = []
    for i in range(row0, row1):  # 行
        for j in range(0, col1):  # 列
            ndvi_list = []
            ndvi_list.append(data1[i][j])
            ndvi_list.append(data2[i][j])
            ndvi_list.append(data3[i][j])
            if get_nan_sum(ndvi_list)>1:
                pass
            else:
                list1.append(ndvi_list)
            ndvi_list = None
    sum_list = get_sum_list(list1)
    list1 = None
    return sum_list


if __name__ == '__main__':
    file_ndvi = r"D:\AAWORK\work\研究方向\研究方向01作物分类内容\风云数据\MERSI-II植被指数旬产品(250M)\NDVI主合成"
    file_out = r"D:\AAWORK\work\2021NDVISUM.tif"
    ndvi_file_list = get_all_file_name(file_ndvi)
    col00, row00, geotrans00, proj00, data_00_ndvi = read_tif(ndvi_file_list[0])
    data_01_ndvi = read_tif02(ndvi_file_list[1])
    data_02_ndvi = read_tif02(ndvi_file_list[2])
    list_row = [0,560,1120,1680,2240,2800,3360,3920,4480,5040,5600]
    for index,i in enumerate(list_row):
        if index <= len(list_row)-2:
            print(list_row[index],list_row[index+1])
            #分段进行操作
            sum_list = get_section(list_row[index],list_row[index+1],col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)
            print(np.array(sum_list))
            # 分段进行保存
            # save_section(sum_list, list_row[index], list_row[index+1], col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)

在这里插入图片描述

3 实现分批保存成TIF文件（所有完整代码）

在2中已经得到了每一批的list结果，我们拿到list结果之后，可以进行保存成tif文件。代码如下：

import os
import numpy as np
from osgeo import gdal, gdalnumeric

def get_sum_list(data_list):
    list1 = []
    for data in data_list:
        sum = 0
        for d in data:
            if not np.isnan(d):
                sum = sum+d
        list1.append(int(sum))
    return list1

def read_tif(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize#图像长度
    row = dataset.RasterYSize#图像宽度
    geotrans = dataset.GetGeoTransform()#读取仿射变换
    proj = dataset.GetProjection()#读取投影
    data = dataset.ReadAsArray()#转为numpy格式
    data = data.astype(np.float32)
    # a = data[0][0]
    # data[data == a] = np.nan
    data[data <= -3] = np.nan
    return [col, row, geotrans, proj, data]

def read_tif02(file):
    data = gdalnumeric.LoadFile(file)
    data = data.astype(np.float32)
    # a = data[0][0]
    # data[data == a] = np.nan
    data[data <= -3] = np.nan
    return data

def get_all_file_name(ndvi_file):
    list1=[]
    if os.path.isdir(ndvi_file):
        fileList = os.listdir(ndvi_file)
        for f in fileList:
            file_name= ndvi_file+"\\"+f
            list1.append(file_name)
        return list1
    else:
        return []

def save_tif(data, file, output):
    """
    保存成tif
    :param data:
    :param file:
    :param output:
    :return:
    """
    ds = gdal.Open(file)
    shape = data.shape
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(output, shape[1], shape[0], 1, gdal.GDT_Float32)#保存的数据类型
    dataset.SetGeoTransform(ds.GetGeoTransform())
    dataset.SetProjection(ds.GetProjection())
    dataset.GetRasterBand(1).WriteArray(data)

def get_nan_sum(ndvi_list):
    """
    得到NAN数据的个数
    :param ndvi_list:
    :return:
    """
    count = 0
    for ndvi in ndvi_list:
        if np.isnan(ndvi):
            count = count+1
    return count

def get_section(row0, row1, col1,data1,data2,data3):
    """
    分段读取数据，读取的数据进行计算
    :param row0:
    :param row1:
    :param col1:
    :param data1:
    :param data2:
    :param data3:
    :return:
    """
    list1 = []
    for i in range(row0, row1):  # 行
        for j in range(0, col1):  # 列
            ndvi_list = []
            ndvi_list.append(data1[i][j])
            ndvi_list.append(data2[i][j])
            ndvi_list.append(data3[i][j])
            if get_nan_sum(ndvi_list)>1:
                pass
            else:
                list1.append(ndvi_list)
            ndvi_list = None
    sum_list = get_sum_list(list1)
    list1 = None
    return sum_list

def save_section(sum_list, row0, row1, col1,data1,data2,data3):
    """
    保存分段的数据
    :param sum_list:
    :param row0:
    :param row1:
    :param col1:
    :param data1:
    :param data2:
    :param data3:
    :return:
    """
    file = r"D:\AAWORK\work\研究方向\研究方向01作物分类内容\风云数据\MERSI-II植被指数旬产品(250M)\kongbai_zhu_250m.tif"#这是一个空白数据，每个像元的值为0
    data = read_tif02(file)
    m = 0
    for i in range(row0, row1):  # 行
        for j in range(0, col1):  # 列
            ndvi_list = []
            ndvi_list.append(data1[i][j])
            ndvi_list.append(data2[i][j])
            ndvi_list.append(data3[i][j])
            if get_nan_sum(ndvi_list)>1:
                pass
            else:
                data[i][j] = sum_list[m]
                m = m + 1
    save_tif(data,file,file_out.replace(".tif","_"+str(row0)+"_"+str(row1)+".tif"))
    data = None


if __name__ == '__main__':
    file_ndvi = r"D:\AAWORK\work\研究方向\研究方向01作物分类内容\风云数据\MERSI-II植被指数旬产品(250M)\NDVI主合成"
    file_out = r"D:\AAWORK\work\2021NDVISUM.tif"
    ndvi_file_list = get_all_file_name(file_ndvi)
    col00, row00, geotrans00, proj00, data_00_ndvi = read_tif(ndvi_file_list[0])
    data_01_ndvi = read_tif02(ndvi_file_list[1])
    data_02_ndvi = read_tif02(ndvi_file_list[2])
    list_row = [0,560,1120,1680,2240,2800,3360,3920,4480,5040,5600]
    for index,i in enumerate(list_row):
        if index <= len(list_row)-2:
            print(list_row[index],list_row[index+1])
            #分段进行操作
            sum_list = get_section(list_row[index],list_row[index+1],col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)
            print(np.array(sum_list))
            # 分段进行保存
            save_section(sum_list, list_row[index], list_row[index+1], col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)

在这里插入图片描述
在这里插入图片描述

4 分段TIF整合到一个TIF

我们要把上述10个TIF文件整合到一个TIF文件里，方法很多，我这里提供一个方法，供大家使用，代码如下：

import os

from osgeo import gdalnumeric, gdal
import numpy as np


def get_all_file_name(file):
    list1=[]
    if os.path.isdir(file):
        fileList = os.listdir(file)
        for f in fileList:
            file_name= file+"\\"+f
            list1.append(file_name)
        return list1
    else:
        return []

def read_tif(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize#图像长度
    row = dataset.RasterYSize#图像宽度
    geotrans = dataset.GetGeoTransform()#读取仿射变换
    proj = dataset.GetProjection()#读取投影
    data = dataset.ReadAsArray()#转为numpy格式
    data = data.astype(np.float32)
    # a = data[0][0]
    # data[data == a] = np.nan
    data[data <= -3] = np.nan
    return [col, row, geotrans, proj, data]

def read_tif02(file):
    data = gdalnumeric.LoadFile(file)
    data = data.astype(np.float32)
    # a = data[0][0]
    # data[data == a] = np.nan
    data[data <= -3] = np.nan
    return data

def save_tif(data, file, output):
    ds = gdal.Open(file)
    shape = data.shape
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(output, shape[1], shape[0], 1, gdal.GDT_Int16)#保存的数据类型
    dataset.SetGeoTransform(ds.GetGeoTransform())
    dataset.SetProjection(ds.GetProjection())
    dataset.GetRasterBand(1).WriteArray(data)

if __name__ == '__main__':
    file_path = r"D:\AAWORK\work\分段数据"
    file_out = r"D:\AAWORK\work\2021NDVISUM.tif"
    file_list = get_all_file_name(file_path)
    data_all = read_tif02(r"D:\AAWORK\work\研究方向\研究方向01作物分类内容\风云数据\MERSI-II植被指数旬产品(250M)\kongbai_zhu_250m.tif")
    for file in file_list:
        data = read_tif02(file)
        data_all = data_all+data
    save_tif(data_all,r"D:\AAWORK\work\研究方向\研究方向01作物分类内容\风云数据\MERSI-II植被指数旬产品(250M)\kongbai_zhu_250m.tif",file_out)

在这里插入图片描述

5 生成一个空白TIF（每个像元值为0的TIF）

思路比较简单，就是遍历每个像元，然后把每个像元的值设置为0，设置为其它可以，然后再进行保存。

from osgeo import gdal
import numpy as np

def read_tif(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize#图像长度
    row = dataset.RasterYSize#图像宽度
    geotrans = dataset.GetGeoTransform()#读取仿射变换
    proj = dataset.GetProjection()#读取投影
    data = dataset.ReadAsArray()#转为numpy格式
    data = data.astype(np.float32)
    # a = data[0][0]
    # data[data == a] = np.nan
    # data[data <= -3] = np.nan
    return [col, row, geotrans, proj, data]

def save_tif(data, file, output):
    ds = gdal.Open(file)
    shape = data.shape
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(output, shape[1], shape[0], 1, gdal.GDT_Int16)#保存的数据类型
    dataset.SetGeoTransform(ds.GetGeoTransform())
    dataset.SetProjection(ds.GetProjection())
    dataset.GetRasterBand(1).WriteArray(data)

if __name__ == '__main__':
    file_path = r"D:\AAWORK\work\2021NDVISUM.tif"
    file_out = r"D:\AAWORK\work\kongbai.tif"
    col, row, geotrans, proj, data = read_tif(file_path)
    for i in range(0,row):
        for j in range(0,col):
            data[i][j] = 0
    
    save_tif(data,file_path,file_out)
本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)

Python图像遥感

Python遥感开发

Python分段读取遥感

Python分段保存遥感

GDAL遥感

Python遥感开发之分段读取和保存遥感数据 的相关文章

随机推荐

  • 编译原理第三版pdf_CS143：编译原理｜PA1：熟悉Cool语言

    本文使用 Zhihu On VSCode 创作并发布 这是本人实现斯坦福CS143变编程作业的笔记 对应第一次作业PA1 环境搭建和一些说明请看 CS143 编译原理 环境搭建HelloWorld 从第一篇搬运一段话 你可能发现了 课程官网

  • 星星之火-56：前传接口 CPRI容器的字长、能力与CPRI速率的对应关系

    CPRI容器的字长 能力与CPRI速率的对应关系 选项 Rate1 Rate2 Rate3 Rate4 Rate5 Rate6 Rate7 Rate8 Rate9 Rate10 Rate倍数 1 2 3 4 5 6 7 8 9 10 码片倍

  • v-if 不生效----vue深入式响应原理

    今天在写的一个行内编辑功能 想通过v if 判断行内是否是编辑状态 来展示是input输入 还是普通显示文本 初始化默认的scope row ipshow1 true 点击编辑按钮 想通过 click scope row ipshow1 f

  • 为什么win11连接wifi频繁掉线？

    如果网络波动比较大的话 就会导致电脑使用过程中不顺畅 网页打开速度都会很缓慢 就有win11用户跟小编反映自己的电脑连接WiFi后总是掉线 非常烦人 这该怎么办 下面就来看看小编为大家整理的几个解决办法 希望可以解决这个问题 方法一 1 首

  • 搭建Ambari Hadoop系统实验

    1 基础环境准备 1 登录到实训系统 检查实训环境 确保有一个CentOS 7系统的虚拟机用来搭建ambari集群 非桌面版虚拟机 首先修改主机名 点击实训页面 虚拟机信息 标签进入非桌面版虚拟机 在 后输入hostnamectl set

  • windows重启nginx服务

    在nginx安装目录下 安全有序停止 nginx exe s quit 启动 start nginx

  • 考研OR工作----计算机操作系统简答题及疑难知识点总结（第三章 处理机调度与死锁）

    上一篇文章总结了一些关于进程的知识点 这章的目的也是根据 计算机操作系统 第四版 汤子瀛 的书来总结一下进程调度和死锁的相关知识点 这一章其实和上一章紧密相连 所以如果没有基础或基础较差 对一些概念还有些模糊 的朋友们先去看上一章的简答题总

  • linux常见面试题

    1 连接linux服务器工具有哪些 SecureCRSecureFX 最好 RealVNC SSHClient putty 比较 SecureCRSecureFX 可以文件传输 可使用命令行 设置字符编码 可开启多个 SSH Client

  • Java拾遗

    这里写自定义目录标题 Hibernate Error Unable to locate persister Caused by java lang NoClassDefFoundError javax sql rowset serial S

  • Nginx超时配置、限流

    目录 一 说明 二 超时配置 三 限流 限制访问频率 限制并发连接数 四 问题记录 五 参考文章 Author Jinwei EditTimes 2020年11月25日17 31 06 一 说明 Nginx 处理的每个请求均有相应的超时设置

  • spring注解及扩展

    1 spring配置注解 spring建议通过注解配置 替代原xml配置方式 使用配置类替代xml配置的优势大体 1 xml配置维护容易出错而且不易检查 java配置类基于java语法检查 对于java程序员更友好 易于维护 2 注解配置

  • NISEDIT如何发布，Qt如何发布文章？难道还有人不会（超详细教学，跟着走，不会你怪我）

    一 自动发布 直接运行即可 不过多阐述 二 手动发布 文件清单 ExamSys exe account txt exam txt Qt5Core dll Qt5Gui dll Qt5Widgets dll libstdc 6 dll lib

  • LeetCode中函数题中“多出来的参数“---returnsize

    转载 关于returnSize 第一次在leetcode上瞎逛就遇到了就遇到了它 int twoSum int nums int numsSize int target int returnSize 1 这个代码的实现并不是什么难解的方法

  • JVM 由哪些部分组成？

    JVM 由哪些部分组成 解析 这是对 JVM 体系结构的考察 答 JVM 的结构基本上由 4 部分组成 类加载器 在 JVM 启动时或者类运行时将需要的 class 加载到 JVM 中 执行引擎 执行引擎的任务是负责执行 class 文件中

  • Zynq7000硬件开发之芯片供电电源功耗（电流）评估

    案头语 单板硬件的主控芯片集成度越来越高 多核处理器越来越多 一块单板可能只需要1块芯片就能满足整体需求 一方面减少设计复杂度 另一面节省PCB面积成本 能同时掌握硬件原理设计以及PCB Layout设计逐渐成为主流 本系列文章同时包含有两

  • ES6详解 快速上手！

    一 Es6 1 1 ES6的概述 ECMAScript的快速发展 编程语言JavaScript是ECMAScript的实现和扩展 ECMAScript是由ECMA 一个类似W3C的标准组织 参与进行标准化的语法规范 ECMAScript定义

  • 【python量化】用python搭建一个股票舆情分析系统

    写在前面 下面的这篇文章将手把手教大家搭建一个简单的股票舆情分析系统 其中将先通过金融界网站爬取指定股票在一段时间的新闻 然后通过百度情感分析接口 用于评估指定股票的正面和反面新闻的占比 以此确定该股票是处于利好还是利空的状态 1 环境准备

  • C++(Liunx) 使用cut截 取出Ubuntu用户的家目录，要求：不能使用“:“作为分割.

    使用cut截 取出Ubuntu用户的家目录 要求 不能使用 作为分割

  • 43.MQ—RabbitMQ

    目录 一 MQ RabbitMQ 1 同步调用与异步调用 1 1 同步调用 1 2 异步调用 2 MQ之间的区别 3 RabbitMQ学习 3 1 docker下载rabbitmq容器 并启动 3 2 RabbitMQ中的几个概念 3 3

  • Python遥感开发之分段读取和保存遥感数据

    Python遥感开发之分段读取和保存遥感数据 1 分段读取数据 2 实现分批读取数据以及进行计算 3 实现分批保存成TIF文件 所有完整代码 4 分段TIF整合到一个TIF 5 生成一个空白TIF 每个像元值为0的TIF 前言 当遇到批量读

热门标签