Python實(shí)現(xiàn)分段讀取和保存遙感數(shù)據(jù)
更新時(shí)間：2023年08月17日 09:56:27 作者：等待著冬天的風(fēng)
當(dāng)遇到批量讀取大量遙感數(shù)據(jù)進(jìn)行運(yùn)算的時(shí)候,如果不進(jìn)行分段讀取操作的話,電腦內(nèi)存可能面臨著不夠使用的情況,所以我們要進(jìn)行分段讀取數(shù)據(jù)然后進(jìn)行運(yùn)算,運(yùn)算結(jié)束之后把這段數(shù)據(jù)保存成tif文件,本文介紹了Python實(shí)現(xiàn)分段讀取和保存遙感數(shù)據(jù),需要的朋友可以參考下
1 分段讀取數(shù)據(jù)
如圖所示，有三個(gè)這樣的數(shù)據(jù)，且該數(shù)據(jù)為5600行6800列，我們可以分成10個(gè)批次分段讀取該TIF數(shù)據(jù)，10個(gè)批次以此為0,560,1120,1680,2240,2800,3360,3920,4480,5040,5600。
代碼實(shí)現(xiàn)：
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()#轉(zhuǎn)為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作物分類內(nèi)容\風(fēng)云數(shù)據(jù)\MERSI-II植被指數(shù)旬產(chǎn)品(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])
#分段進(jìn)行操作
# sum_list = get_section(list_row[index],list_row[index+1],col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)
# 分段進(jìn)行保存
# save_section(sum_list, list_row[index], list_row[index+1], col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)
2 實(shí)現(xiàn)分批讀取數(shù)據(jù)以及進(jìn)行計(jì)算
拿到開始的行和結(jié)束的行數(shù)，進(jìn)行分批讀取數(shù)據(jù)并進(jìn)行計(jì)算，（這里求和求的是整數(shù)，如有需要的話可以自己更改）代碼如下：
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()#轉(zhuǎn)為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數(shù)據(jù)的個(gè)數(shù)
: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):
"""
分段讀取數(shù)據(jù)，讀取的數(shù)據(jù)進(jìn)行計(jì)算
: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作物分類內(nèi)容\風(fēng)云數(shù)據(jù)\MERSI-II植被指數(shù)旬產(chǎn)品(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])
#分段進(jìn)行操作
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))
# 分段進(jìn)行保存
# save_section(sum_list, list_row[index], list_row[index+1], col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)
3 實(shí)現(xiàn)分批保存成TIF文件（所有完整代碼）
在2中已經(jīng)得到了每一批的list結(jié)果，我們拿到list結(jié)果之后，可以進(jìn)行保存成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()#轉(zhuǎn)為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)#保存的數(shù)據(jù)類型
dataset.SetGeoTransform(ds.GetGeoTransform())
dataset.SetProjection(ds.GetProjection())
dataset.GetRasterBand(1).WriteArray(data)
def get_nan_sum(ndvi_list):
"""
得到NAN數(shù)據(jù)的個(gè)數(shù)
: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):
"""
分段讀取數(shù)據(jù)，讀取的數(shù)據(jù)進(jìn)行計(jì)算
: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):
"""
保存分段的數(shù)據(jù)
:param sum_list:
:param row0:
:param row1:
:param col1:
:param data1:
:param data2:
:param data3:
:return:
"""
file = r"D:\AAWORK\work\研究方向\研究方向01作物分類內(nèi)容\風(fēng)云數(shù)據(jù)\MERSI-II植被指數(shù)旬產(chǎn)品(250M)\kongbai_zhu_250m.tif"#這是一個(gè)空白數(shù)據(jù)，每個(gè)像元的值為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作物分類內(nèi)容\風(fēng)云數(shù)據(jù)\MERSI-II植被指數(shù)旬產(chǎn)品(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])
#分段進(jìn)行操作
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))
# 分段進(jìn)行保存
save_section(sum_list, list_row[index], list_row[index+1], col00,data_00_ndvi,data_01_ndvi,data_02_ndvi)
4 分段TIF整合到一個(gè)TIF
我們要把上述10個(gè)TIF文件整合到一個(gè)TIF文件里，方法很多，我這里提供一個(gè)方法，供大家使用，代碼如下：
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()#轉(zhuǎn)為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)#保存的數(shù)據(jù)類型
dataset.SetGeoTransform(ds.GetGeoTransform())
dataset.SetProjection(ds.GetProjection())
dataset.GetRasterBand(1).WriteArray(data)
if __name__ == '__main__':
file_path = r"D:\AAWORK\work\分段數(shù)據(jù)"
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作物分類內(nèi)容\風(fēng)云數(shù)據(jù)\MERSI-II植被指數(shù)旬產(chǎn)品(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作物分類內(nèi)容\風(fēng)云數(shù)據(jù)\MERSI-II植被指數(shù)旬產(chǎn)品(250M)\kongbai_zhu_250m.tif",file_out)
5 生成一個(gè)空白TIF（每個(gè)像元值為0的TIF）
思路比較簡單，就是遍歷每個(gè)像元，然后把每個(gè)像元的值設(shè)置為0，設(shè)置為其它可以，然后再進(jìn)行保存。
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()#轉(zhuǎn)為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)#保存的數(shù)據(jù)類型
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)以上就是Python實(shí)現(xiàn)分段讀取和保存遙感數(shù)據(jù)的詳細(xì)內(nèi)容，更多關(guān)于Python遙感數(shù)據(jù)的資料請(qǐng)關(guān)注腳本之家其它相關(guān)文章！
相關(guān)文章
這篇文章主要介紹了Python中使用裝飾器時(shí)需要注意的一些問題,裝飾器是Python學(xué)習(xí)進(jìn)階中的重要知識(shí)點(diǎn),需要的朋友可以參考下2015-05-05
在本文里小編給讀者們分享一篇關(guān)于pycharm怎么新建一個(gè)python工程的知識(shí)點(diǎn)和步驟內(nèi)容，需要的朋友們學(xué)習(xí)下。2019-07-07
win10環(huán)境下python3.5安裝步驟圖文教程
本文通過圖文并茂的形式給大家介紹了win10環(huán)境下python3.5安裝步驟,需要的朋友可以參考下2017-02-02
使用django-guardian實(shí)現(xiàn)django-admin的行級(jí)權(quán)限控制的方法
這篇文章主要介紹了使用django-guardian實(shí)現(xiàn)django-admin的行級(jí)權(quán)限控制的方法，小編覺得挺不錯(cuò)的，現(xiàn)在分享給大家，也給大家做個(gè)參考。一起跟隨小編過來看看吧2018-10-10
Python?DataFrame處理缺失值的完整指南與實(shí)戰(zhàn)技巧
在數(shù)據(jù)分析工作中,缺失值（Missing?Values）是不可避免的挑戰(zhàn),本文將系統(tǒng)介紹DataFrame缺失值的識(shí)別、處理策略和實(shí)戰(zhàn)技巧,希望對(duì)大家有所幫助2026-02-02
Python實(shí)現(xiàn)partial改變方法默認(rèn)參數(shù)
這篇文章主要介紹了Python實(shí)現(xiàn)partial改變方法默認(rèn)參數(shù),需要的朋友可以參考下2014-08-08
re模塊是Python中使用正則表達(dá)式的最基礎(chǔ)的模塊,re模塊的這些功能覆蓋了正則表達(dá)式的常見用法,使用re模塊可以簡化字符串的模式匹配、信息提取、過濾替換、切分等操作,本文給大家介紹正則表達(dá)式模塊常用方法,感興趣的朋友跟隨小編一起看看吧2023-08-08
Python?Flask中Cookie和Session區(qū)別詳解
Flask是一個(gè)使用?Python?編寫的輕量級(jí)?Web?應(yīng)用框架。其?WSGI?工具箱采用?Werkzeug?，模板引擎則使用?Jinja2?。Flask使用?BSD?授權(quán)。Flask也被稱為?“microframework”?，因?yàn)樗褂煤唵蔚暮诵?，?extension?增加其他功能，F(xiàn)lask中Cookie和Session有什么區(qū)別呢2022-07-07
python3 requests庫實(shí)現(xiàn)多圖片爬取教程
今天小編就為大家分享一篇python3 requests庫實(shí)現(xiàn)多圖片爬取教程，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過來看看吧2019-12-12
Blender Python編程實(shí)現(xiàn)程序化建模生成超形示例詳解
這篇文章主要為大家介紹了Blender Python編程實(shí)現(xiàn)程序化建模生成超形示例詳解，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進(jìn)步，早日升職加薪2022-08-08