Zunächst müssen Sie die verwendeten Module und verschiedene Pfade zum Speichern von Rasterbildern vorbereiten.
import os import copy import numpy as np import pylab as plt from osgeo import gdal # rt_file_path="G:/Postgraduate/LAI_Glass_RTlab/Rc_Lai_A2018161_h22v03.tif" # gl_file_path="G:/Postgraduate/LAI_Glass_RTlab/GLASS01E01.V50.A2018161.h22v03.2020323.hdf" # out_file_path="G:/Postgraduate/LAI_Glass_RTlab/test.tif" rt_file_path="I:/LAI_RTLab/A2018161/" gl_file_path="I:/LAI_Glass/2018161/" out_file_path="I:/LAI_Dif/"
Unter diesen ist rt_file_path der Speicherpfad seiner eigenen Produkte, gl_file_path ist der Speicherpfad von GLASS Produkte, < code>out_file_path ist der Speicherpfad des Endergebnisses nach der Differenzverarbeitung zwischen den beiden Rastern. rt_file_path为自有产品的存放路径，gl_file_path为GLASS产品的存放路径，out_file_path为最终二者栅格做完差值处理后结果的存放路径。
接下来，需要将全部待处理的栅格图像用os.listdir()进行获取，并用for循环进行循环批量处理操作的准备。
rt_file_list=os.listdir(rt_file_path)
for rt_file in rt_file_list:
file_name_split=rt_file.split("_")
rt_hv=file_name_split[3][:-4]
gl_file_list=os.listdir(gl_file_path)
for gl_file in gl_file_list:
if rt_hv in gl_file:
rt_file_tif_path=rt_file_path+rt_file
gl_file_tif_path=gl_file_path+gl_file其中，由于本文需求是对两种产品做差，因此首先需要结合二者的hv分幅编号，将同一分幅编号的两景遥感影像放在一起；因此，依据自有产品文件名的特征，选择.split()进行字符串分割，并随后截取获得遥感影像的hv分幅编号。
前述1.1部分已经配置好了输出文件存放的路径，但是还没有进行输出文件文件名的配置；因此这里我们需要配置好每一个做差后的遥感影像的文件存放路径与名称。其中，我们就直接以遥感影像的hv编号作为输出结果文件名。
DRT_out_file_path=out_file_path+"DRT/"
if not os.path.exists(DRT_out_file_path):
os.makedirs(DRT_out_file_path)
DRT_out_file_tif_path=os.path.join(DRT_out_file_path,rt_hv+".tif")
eco_out_file_path=out_file_path+"eco/"
if not os.path.exists(eco_out_file_path):
os.makedirs(eco_out_file_path)
eco_out_file_tif_path=os.path.join(eco_out_file_path,rt_hv+".tif")
wat_out_file_path=out_file_path+"wat/"
if not os.path.exists(wat_out_file_path):
os.makedirs(wat_out_file_path)
wat_out_file_tif_path=os.path.join(wat_out_file_path,rt_hv+".tif")
tim_out_file_path=out_file_path+"tim/"
if not os.path.exists(tim_out_file_path):
os.makedirs(tim_out_file_path)
tim_out_file_tif_path=os.path.join(tim_out_file_path,rt_hv+".tif")这一部分代码分为了四个部分，是因为自有产品的LAI是分别依据四种算法得到的，在做差时需要每一种算法分别和GLASS产品进行相减，因此配置了四个输出路径文件夹。
接下来，利用gdal模块对.tif与.hdf等两种栅格图像加以读取。
rt_raster=gdal.Open(rt_file_path+rt_file)
rt_band_num=rt_raster.RasterCount
rt_raster_array=rt_raster.ReadAsArray()
rt_lai_array=rt_raster_array[0]
rt_qa_array=rt_raster_array[1]
rt_lai_band=rt_raster.GetRasterBand(1)
# rt_lai_nodata=rt_lai_band.GetNoDataValue()
# rt_lai_nodata=32767
# rt_lai_mask=np.ma.masked_equal(rt_lai_array,rt_lai_nodata)
rt_lai_array_mask=np.where(rt_lai_array>30000,np.nan,rt_lai_array)
rt_lai_array_fin=rt_lai_array_mask*0.001
gl_raster=gdal.Open(gl_file_path+gl_file)
gl_band_num=gl_raster.RasterCount
gl_raster_array=gl_raster.ReadAsArray()
gl_lai_array=gl_raster_array
gl_lai_band=gl_raster.GetRasterBand(1)
gl_lai_array_mask=np.where(gl_lai_array>1000,np.nan,gl_lai_array)
gl_lai_array_fin=gl_lai_array_mask*0.01
row=rt_raster.RasterYSize
col=rt_raster.RasterXSize
geotransform=rt_raster.GetGeoTransform()
projection=rt_raster.GetProjection()首先，以上述代码的第一段为例进行讲解。其中，gdal.Open()读取栅格图像；.RasterCount获取栅格图像波段数量；.ReadAsArray()将栅格图像各波段的信息读取为Array格式，当波段数量大于1时，其共有三维，第一维为波段的个数；rt_raster_array[0]表示取Array中的第一个波段，在本文中也就是自有产品的LAI波段；rt_qa_array=rt_raster_array[1]则表示取出第二个波段，在本文中也就是自有产品的QA波段；.GetRasterBand(1)表示获取栅格图像中的第一个波段（注意，这里序号不是从0开始而是从1开始）；np.where(rt_lai_array>30000,np.nan,rt_lai_array)表示利用np.where()函数对Array中第一个波段中像素>30000加以选取，并将其设置为nan，其他值不变。这一步骤是消除图像中填充值、Nodata值的方法。最后一句*0.001是将图层原有的缩放系数复原。
其次，上述代码第三段为获取栅格行、列数与投影变换信息。
接下来，首先对自有产品与GLASS产品加以做差操作，随后需要对四种算法分别加以提取。
lai_dif=rt_lai_array_fin-gl_lai_array_fin
lai_dif=lai_dif*1000
rt_qa_array_bin=copy.copy(rt_qa_array)
rt_qa_array_row,rt_qa_array_col=rt_qa_array.shape
for i in range(rt_qa_array_row):
for j in range(rt_qa_array_col):
rt_qa_array_bin[i][j]="{:012b}".format(rt_qa_array_bin[i][j])[-4:]
# DRT_pixel_pos=np.where((rt_qa_array_bin>=100) & (rt_qa_array_bin==11))
# eco_pixel_pos=np.where((rt_qa_array_bin<100) & (rt_qa_array_bin==111))
# wat_pixel_pos=np.where((rt_qa_array_bin<1000) & (rt_qa_array_bin==1011))
# tim_pixel_pos=np.where((rt_qa_array_bin<1100) & (rt_qa_array_bin==1111))
# colormap=plt.cm.Greens
# plt.figure(1)
# # plt.subplot(2,4,1)
# plt.imshow(rt_lai_array_fin,cmap=colormap,interpolation='none')
# plt.title("RT_LAI")
# plt.colorbar()
# plt.figure(2)
# # plt.subplot(2,4,2)
# plt.imshow(gl_lai_array_fin,cmap=colormap,interpolation='none')
# plt.title("GLASS_LAI")
# plt.colorbar()
# plt.figure(3)
# dif_colormap=plt.cm.get_cmap("Spectral")
# plt.imshow(lai_dif,cmap=dif_colormap,interpolation='none')
# plt.title("Difference_LAI (RT-GLASS)")
# plt.colorbar()
DRT_lai_dif_array=np.where((rt_qa_array_bin>=100) | (rt_qa_array_bin==11),
np.nan,lai_dif)
eco_lai_dif_array=np.where((rt_qa_array_bin<100) | (rt_qa_array_bin==111),
np.nan,lai_dif)
wat_lai_dif_array=np.where((rt_qa_array_bin<1000) | (rt_qa_array_bin==1011),
np.nan,lai_dif)
tim_lai_dif_array=np.where((rt_qa_array_bin<1100) | (rt_qa_array_bin==1111),
np.nan,lai_dif)
# plt.figure(4)
# plt.imshow(DRT_lai_dif_array)
# plt.colorbar()
# plt.figure(5)
# plt.imshow(eco_lai_dif_array)
# plt.colorbar()
# plt.figure(6)
# plt.imshow(wat_lai_dif_array)
# plt.colorbar()
# plt.figure(7)
# plt.imshow(tim_lai_dif_array)
# plt.colorbar()其中，上述代码前两句为差值计算与数据化整。将数据转换为整数，可以减少结果数据图层的数据量（因为不需要存储小数了）。
随后，开始依据QA波段进行数据筛选与掩膜。其实各类遥感影像（例如MODIS、Landsat等）的QA波段都是比较近似的：通过一串二进制码来表示遥感影像的质量、信息等，其中不同的比特位往往都代表着一种特性。例如下图所示为Landsat Collection 2 Level-2的QA波段含义。
在这里，QA波段原本为十进制（一般遥感影像为了节省空间，QA波段都是写成十进制的形式），因此需要将其转换为二进制；随后通过获取指定需要的二进制数据位数（在本文中也就是能确定自有产品中这一像素来自于哪一种算法的二进制位数），从而判断这一像素所得LAI是通过哪一种算法得到的，从而将每种算法对应的像素分别放在一起处理。DRT_lai_dif_array等四个变量分别表示四种算法中，除了自己这一种算法得到的像素之外的其他所有像素；之所以选择这种方式，是因为后期我们可以将其直接掩膜掉，那么剩下的就是这种算法自身的像素了。
其中，上述代码注释掉的plt
os.listdir()</code >Acquire and use <code>for< verwenden /code>-Schleife zur Vorbereitung der Schleifenstapelverarbeitung. #🎜🎜#<div class="code" style="position:relative; padding:0px; margin:0px;"><div class="code" style="position:relative; padding:0px; margin:0px;"><pre class="brush:py;"> driver=gdal.GetDriverByName("Gtiff")
out_DRT_lai=driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_DRT_lai.SetGeoTransform(geotransform)
out_DRT_lai.SetProjection(projection)
out_DRT_lai.GetRasterBand(1).WriteArray(DRT_lai_dif_array)
out_DRT_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_eco_lai.SetGeoTransform(geotransform)
out_eco_lai.SetProjection(projection)
out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
out_eco_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_wat_lai.SetGeoTransform(geotransform)
out_wat_lai.SetProjection(projection)
out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
out_wat_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_tim_lai.SetGeoTransform(geotransform)
out_tim_lai.SetProjection(projection)
out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
out_tim_lai=None
print(rt_hv)</pre><div class="contentsignin">Nach dem Login kopieren</div></div><div class="contentsignin">Nach dem Login kopieren</div></div>#🎜🎜# Da die Anforderung dieses Artikels darin besteht, die beiden Produkte zu vergleichen, müssen zunächst die <code>hv-Rahmennummern der beiden kombiniert werden, um die beiden Fernbedienungen zu kombinieren Wählen Sie daher basierend auf den Eigenschaften des eigenen Produktdateinamens .split() aus, um eine Zeichenfolgenaufteilung durchzuführen, und fangen Sie sie dann ab, um die Bilder zu erfassen hv-Frame der Seriennummer des Fernerkundungsbildes. #🎜🎜##🎜🎜#1.3 Vorbereitung des Ausgabedateinamens #🎜🎜##🎜🎜#Der Pfad zum Speichern der Ausgabedatei wurde im oben genannten Abschnitt 1.1 konfiguriert, der Ausgabedateiname jedoch schon Die Konfiguration ist noch nicht erfolgt. Daher müssen wir hier den Dateispeicherpfad und den Namen jedes Fernerkundungsbildes nach der Verarbeitung konfigurieren. Unter anderem verwenden wir direkt die hv-Nummer des Fernerkundungsbildes als Namen der Ausgabeergebnisdatei. #🎜🎜## -*- coding: utf-8 -*-
"""
Created on Thu Jul 15 19:36:15 2021
@author: fkxxgis
"""
import os
import copy
import numpy as np
import pylab as plt
from osgeo import gdal
# rt_file_path="G:/Postgraduate/LAI_Glass_RTlab/Rc_Lai_A2018161_h22v03.tif"
# gl_file_path="G:/Postgraduate/LAI_Glass_RTlab/GLASS01E01.V50.A2018161.h22v03.2020323.hdf"
# out_file_path="G:/Postgraduate/LAI_Glass_RTlab/test.tif"
rt_file_path="I:/LAI_RTLab/A2018161/"
gl_file_path="I:/LAI_Glass/2018161/"
out_file_path="I:/LAI_Dif/"
rt_file_list=os.listdir(rt_file_path)
for rt_file in rt_file_list:
file_name_split=rt_file.split("_")
rt_hv=file_name_split[3][:-4]
gl_file_list=os.listdir(gl_file_path)
for gl_file in gl_file_list:
if rt_hv in gl_file:
rt_file_tif_path=rt_file_path+rt_file
gl_file_tif_path=gl_file_path+gl_file
DRT_out_file_path=out_file_path+"DRT/"
if not os.path.exists(DRT_out_file_path):
os.makedirs(DRT_out_file_path)
DRT_out_file_tif_path=os.path.join(DRT_out_file_path,rt_hv+".tif")
eco_out_file_path=out_file_path+"eco/"
if not os.path.exists(eco_out_file_path):
os.makedirs(eco_out_file_path)
eco_out_file_tif_path=os.path.join(eco_out_file_path,rt_hv+".tif")
wat_out_file_path=out_file_path+"wat/"
if not os.path.exists(wat_out_file_path):
os.makedirs(wat_out_file_path)
wat_out_file_tif_path=os.path.join(wat_out_file_path,rt_hv+".tif")
tim_out_file_path=out_file_path+"tim/"
if not os.path.exists(tim_out_file_path):
os.makedirs(tim_out_file_path)
tim_out_file_tif_path=os.path.join(tim_out_file_path,rt_hv+".tif")
rt_raster=gdal.Open(rt_file_path+rt_file)
rt_band_num=rt_raster.RasterCount
rt_raster_array=rt_raster.ReadAsArray()
rt_lai_array=rt_raster_array[0]
rt_qa_array=rt_raster_array[1]
rt_lai_band=rt_raster.GetRasterBand(1)
# rt_lai_nodata=rt_lai_band.GetNoDataValue()
# rt_lai_nodata=32767
# rt_lai_mask=np.ma.masked_equal(rt_lai_array,rt_lai_nodata)
rt_lai_array_mask=np.where(rt_lai_array>30000,np.nan,rt_lai_array)
rt_lai_array_fin=rt_lai_array_mask*0.001
gl_raster=gdal.Open(gl_file_path+gl_file)
gl_band_num=gl_raster.RasterCount
gl_raster_array=gl_raster.ReadAsArray()
gl_lai_array=gl_raster_array
gl_lai_band=gl_raster.GetRasterBand(1)
gl_lai_array_mask=np.where(gl_lai_array>1000,np.nan,gl_lai_array)
gl_lai_array_fin=gl_lai_array_mask*0.01
row=rt_raster.RasterYSize
col=rt_raster.RasterXSize
geotransform=rt_raster.GetGeoTransform()
projection=rt_raster.GetProjection()
lai_dif=rt_lai_array_fin-gl_lai_array_fin
lai_dif=lai_dif*1000
rt_qa_array_bin=copy.copy(rt_qa_array)
rt_qa_array_row,rt_qa_array_col=rt_qa_array.shape
for i in range(rt_qa_array_row):
for j in range(rt_qa_array_col):
rt_qa_array_bin[i][j]="{:012b}".format(rt_qa_array_bin[i][j])[-4:]
# DRT_pixel_pos=np.where((rt_qa_array_bin>=100) & (rt_qa_array_bin==11))
# eco_pixel_pos=np.where((rt_qa_array_bin<100) & (rt_qa_array_bin==111))
# wat_pixel_pos=np.where((rt_qa_array_bin<1000) & (rt_qa_array_bin==1011))
# tim_pixel_pos=np.where((rt_qa_array_bin<1100) & (rt_qa_array_bin==1111))
# colormap=plt.cm.Greens
# plt.figure(1)
# # plt.subplot(2,4,1)
# plt.imshow(rt_lai_array_fin,cmap=colormap,interpolation='none')
# plt.title("RT_LAI")
# plt.colorbar()
# plt.figure(2)
# # plt.subplot(2,4,2)
# plt.imshow(gl_lai_array_fin,cmap=colormap,interpolation='none')
# plt.title("GLASS_LAI")
# plt.colorbar()
# plt.figure(3)
# dif_colormap=plt.cm.get_cmap("Spectral")
# plt.imshow(lai_dif,cmap=dif_colormap,interpolation='none')
# plt.title("Difference_LAI (RT-GLASS)")
# plt.colorbar()
DRT_lai_dif_array=np.where((rt_qa_array_bin>=100) | (rt_qa_array_bin==11),
np.nan,lai_dif)
eco_lai_dif_array=np.where((rt_qa_array_bin<100) | (rt_qa_array_bin==111),
np.nan,lai_dif)
wat_lai_dif_array=np.where((rt_qa_array_bin<1000) | (rt_qa_array_bin==1011),
np.nan,lai_dif)
tim_lai_dif_array=np.where((rt_qa_array_bin<1100) | (rt_qa_array_bin==1111),
np.nan,lai_dif)
# plt.figure(4)
# plt.imshow(DRT_lai_dif_array)
# plt.colorbar()
# plt.figure(5)
# plt.imshow(eco_lai_dif_array)
# plt.colorbar()
# plt.figure(6)
# plt.imshow(wat_lai_dif_array)
# plt.colorbar()
# plt.figure(7)
# plt.imshow(tim_lai_dif_array)
# plt.colorbar()
driver=gdal.GetDriverByName("Gtiff")
out_DRT_lai=driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_DRT_lai.SetGeoTransform(geotransform)
out_DRT_lai.SetProjection(projection)
out_DRT_lai.GetRasterBand(1).WriteArray(DRT_lai_dif_array)
out_DRT_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_eco_lai.SetGeoTransform(geotransform)
out_eco_lai.SetProjection(projection)
out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
out_eco_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_wat_lai.SetGeoTransform(geotransform)
out_wat_lai.SetProjection(projection)
out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
out_wat_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_tim_lai.SetGeoTransform(geotransform)
out_tim_lai.SetProjection(projection)
out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
out_tim_lai=None
print(rt_hv)gdal, um .tif und < Code zu vergleichen >.hdf und zwei weitere Rasterbilder können gelesen werden. #🎜🎜#rrreee#🎜🎜#Nehmen wir zunächst den ersten Absatz des obigen Codes als Beispiel. Unter anderem liest gdal.Open() das Rasterbild; .RasterCount ruft die Anzahl der Rasterbildbänder ab; .ReadAsArray() konvertiert das Rasterbild in Die Informationen jedes Bandes des Bildes werden im Array-Format gelesen. Wenn die Anzahl der Bänder größer als 1 ist, hat es drei Dimensionen und die erste Dimension ist die Anzahl der Bänder; rt_raster_array[ 0] bedeutet, dass wir das erste Band im Array nehmen, das in diesem Artikel unser eigenes LAI-Band ist product; rt_qa_array=rt_raster_array[1 ] bedeutet, das zweite Band herauszunehmen, das in diesem Artikel das QA-Band unseres eigenen Produkts .GetRasterBand(1) ist ) bedeutet, das Rasterbild zu erhalten. Das erste Band (beachten Sie, dass die Seriennummer hier nicht bei 0, sondern bei 1 beginnt); where(rt_lai_array>30000,np.nan, rt_lai_array) bedeutet, dass die Funktion np.where() verwendet wird, um das Pixel >30000 im ersten Band auszuwählen in Array und setzen Sie es auf nan, wobei andere Werte unverändert bleiben. Dieser Schritt besteht darin, Auffüllungen und Nodata-Werte aus dem Bild zu entfernen. Der letzte Satz *0.001 dient dazu, den ursprünglichen Skalierungsfaktor der Ebene wiederherzustellen. #🎜🎜##🎜🎜#Zweitens besteht der dritte Abschnitt des obigen Codes darin, die Rasterzeilen-, Spaltennummer- und Projektionstransformationsinformationen abzurufen. #🎜🎜##🎜🎜#1.5 Differenzberechnung und QA-Band-Screening #🎜🎜##🎜🎜# Als nächstes führen wir zunächst eine Differenzoperation an unseren eigenen Produkten und GLAS-Produkten durch, und dann müssen wir dies tun Vier Algorithmen werden separat extrahiert. #🎜🎜#rrreee#🎜🎜# Darunter sind die ersten beiden Sätze des obigen Codes Differenzberechnung und Datenrundung. Durch die Konvertierung der Daten in Ganzzahlen wird die Datenmenge in der resultierenden Datenschicht reduziert (da keine Dezimalzahlen gespeichert werden müssen). #🎜🎜##🎜🎜#Dann begann die Datenüberprüfung und -maskierung basierend auf dem QA-Band. Tatsächlich sind die QA-Bänder verschiedener Fernerkundungsbilder (wie MODIS, Landsat usw.) relativ ähnlich: Fernerkundung wird dargestellt durch eine Folge binärer Codes für Bildqualität, Informationen usw., wobei verschiedene Bits häufig ein Merkmal darstellen. Die folgende Abbildung zeigt beispielsweise die Bedeutung des QA-Bands von Landsat Collection 2 Level-2. #🎜🎜##🎜🎜#Hier war das QA-Band ursprünglich in Dezimalform (um in allgemeinen Fernerkundungsbildern Platz zu sparen, ist das QA-Band eingeschrieben Dezimalform), also ist es notwendig, es in eine Binärform umzuwandeln; dann wird die Anzahl der erforderlichen Binärdatenbits ermittelt (in diesem Artikel ist es die Anzahl der Binärbits, die bestimmen kann, in welchem Algorithmus dieses Pixel verwendet wird). Ihr eigenes Produkt stammt). Dies bestimmt, durch welchen Algorithmus der für dieses Pixel erhaltene LAI erhalten wurde, und die Pixel, die jedem Algorithmus entsprechen, werden zusammen verarbeitet. Vier Variablen, darunter DRT_lai_dif_array, stellen jeweils alle Pixel in den vier Algorithmen dar, mit Ausnahme der vom eigenen Algorithmus erhaltenen Pixel. Der Grund für die Wahl dieser Methode besteht darin, dass wir sie später direkt verwenden können. Übrig bleiben die Pixel dieses Algorithmus selbst. #🎜🎜##🎜🎜# Unter anderem kann der im obigen Code auskommentierte plt-Inhalt zum Zeichnen räumlicher Verteilungskarten verwendet werden. Wenn Sie interessiert sind, können Sie versuchen, ihn zu verwenden. #🎜🎜#接下来，将我们完成上述差值计算与依据算法进行筛选后的图像保存。
driver=gdal.GetDriverByName("Gtiff")
out_DRT_lai=driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_DRT_lai.SetGeoTransform(geotransform)
out_DRT_lai.SetProjection(projection)
out_DRT_lai.GetRasterBand(1).WriteArray(DRT_lai_dif_array)
out_DRT_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_eco_lai.SetGeoTransform(geotransform)
out_eco_lai.SetProjection(projection)
out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
out_eco_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_wat_lai.SetGeoTransform(geotransform)
out_wat_lai.SetProjection(projection)
out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
out_wat_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_tim_lai.SetGeoTransform(geotransform)
out_tim_lai.SetProjection(projection)
out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
out_tim_lai=None
print(rt_hv)其中，.GetDriverByName("Gtiff")表示保存为.tif格式的GeoTIFF文件；driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)表示按照路径、行列数、波段数与数据格式等建立一个新的栅格图层，作为输出图层的框架；其后表示分别将地理投影转换信息与像素具体数值分别赋予这一新建的栅格图层；最后=None表示将其从内存空间中释放，完成写入与保存工作。
本文所需完整代码如下：
# -*- coding: utf-8 -*-
"""
Created on Thu Jul 15 19:36:15 2021
@author: fkxxgis
"""
import os
import copy
import numpy as np
import pylab as plt
from osgeo import gdal
# rt_file_path="G:/Postgraduate/LAI_Glass_RTlab/Rc_Lai_A2018161_h22v03.tif"
# gl_file_path="G:/Postgraduate/LAI_Glass_RTlab/GLASS01E01.V50.A2018161.h22v03.2020323.hdf"
# out_file_path="G:/Postgraduate/LAI_Glass_RTlab/test.tif"
rt_file_path="I:/LAI_RTLab/A2018161/"
gl_file_path="I:/LAI_Glass/2018161/"
out_file_path="I:/LAI_Dif/"
rt_file_list=os.listdir(rt_file_path)
for rt_file in rt_file_list:
file_name_split=rt_file.split("_")
rt_hv=file_name_split[3][:-4]
gl_file_list=os.listdir(gl_file_path)
for gl_file in gl_file_list:
if rt_hv in gl_file:
rt_file_tif_path=rt_file_path+rt_file
gl_file_tif_path=gl_file_path+gl_file
DRT_out_file_path=out_file_path+"DRT/"
if not os.path.exists(DRT_out_file_path):
os.makedirs(DRT_out_file_path)
DRT_out_file_tif_path=os.path.join(DRT_out_file_path,rt_hv+".tif")
eco_out_file_path=out_file_path+"eco/"
if not os.path.exists(eco_out_file_path):
os.makedirs(eco_out_file_path)
eco_out_file_tif_path=os.path.join(eco_out_file_path,rt_hv+".tif")
wat_out_file_path=out_file_path+"wat/"
if not os.path.exists(wat_out_file_path):
os.makedirs(wat_out_file_path)
wat_out_file_tif_path=os.path.join(wat_out_file_path,rt_hv+".tif")
tim_out_file_path=out_file_path+"tim/"
if not os.path.exists(tim_out_file_path):
os.makedirs(tim_out_file_path)
tim_out_file_tif_path=os.path.join(tim_out_file_path,rt_hv+".tif")
rt_raster=gdal.Open(rt_file_path+rt_file)
rt_band_num=rt_raster.RasterCount
rt_raster_array=rt_raster.ReadAsArray()
rt_lai_array=rt_raster_array[0]
rt_qa_array=rt_raster_array[1]
rt_lai_band=rt_raster.GetRasterBand(1)
# rt_lai_nodata=rt_lai_band.GetNoDataValue()
# rt_lai_nodata=32767
# rt_lai_mask=np.ma.masked_equal(rt_lai_array,rt_lai_nodata)
rt_lai_array_mask=np.where(rt_lai_array>30000,np.nan,rt_lai_array)
rt_lai_array_fin=rt_lai_array_mask*0.001
gl_raster=gdal.Open(gl_file_path+gl_file)
gl_band_num=gl_raster.RasterCount
gl_raster_array=gl_raster.ReadAsArray()
gl_lai_array=gl_raster_array
gl_lai_band=gl_raster.GetRasterBand(1)
gl_lai_array_mask=np.where(gl_lai_array>1000,np.nan,gl_lai_array)
gl_lai_array_fin=gl_lai_array_mask*0.01
row=rt_raster.RasterYSize
col=rt_raster.RasterXSize
geotransform=rt_raster.GetGeoTransform()
projection=rt_raster.GetProjection()
lai_dif=rt_lai_array_fin-gl_lai_array_fin
lai_dif=lai_dif*1000
rt_qa_array_bin=copy.copy(rt_qa_array)
rt_qa_array_row,rt_qa_array_col=rt_qa_array.shape
for i in range(rt_qa_array_row):
for j in range(rt_qa_array_col):
rt_qa_array_bin[i][j]="{:012b}".format(rt_qa_array_bin[i][j])[-4:]
# DRT_pixel_pos=np.where((rt_qa_array_bin>=100) & (rt_qa_array_bin==11))
# eco_pixel_pos=np.where((rt_qa_array_bin<100) & (rt_qa_array_bin==111))
# wat_pixel_pos=np.where((rt_qa_array_bin<1000) & (rt_qa_array_bin==1011))
# tim_pixel_pos=np.where((rt_qa_array_bin<1100) & (rt_qa_array_bin==1111))
# colormap=plt.cm.Greens
# plt.figure(1)
# # plt.subplot(2,4,1)
# plt.imshow(rt_lai_array_fin,cmap=colormap,interpolation='none')
# plt.title("RT_LAI")
# plt.colorbar()
# plt.figure(2)
# # plt.subplot(2,4,2)
# plt.imshow(gl_lai_array_fin,cmap=colormap,interpolation='none')
# plt.title("GLASS_LAI")
# plt.colorbar()
# plt.figure(3)
# dif_colormap=plt.cm.get_cmap("Spectral")
# plt.imshow(lai_dif,cmap=dif_colormap,interpolation='none')
# plt.title("Difference_LAI (RT-GLASS)")
# plt.colorbar()
DRT_lai_dif_array=np.where((rt_qa_array_bin>=100) | (rt_qa_array_bin==11),
np.nan,lai_dif)
eco_lai_dif_array=np.where((rt_qa_array_bin<100) | (rt_qa_array_bin==111),
np.nan,lai_dif)
wat_lai_dif_array=np.where((rt_qa_array_bin<1000) | (rt_qa_array_bin==1011),
np.nan,lai_dif)
tim_lai_dif_array=np.where((rt_qa_array_bin<1100) | (rt_qa_array_bin==1111),
np.nan,lai_dif)
# plt.figure(4)
# plt.imshow(DRT_lai_dif_array)
# plt.colorbar()
# plt.figure(5)
# plt.imshow(eco_lai_dif_array)
# plt.colorbar()
# plt.figure(6)
# plt.imshow(wat_lai_dif_array)
# plt.colorbar()
# plt.figure(7)
# plt.imshow(tim_lai_dif_array)
# plt.colorbar()
driver=gdal.GetDriverByName("Gtiff")
out_DRT_lai=driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_DRT_lai.SetGeoTransform(geotransform)
out_DRT_lai.SetProjection(projection)
out_DRT_lai.GetRasterBand(1).WriteArray(DRT_lai_dif_array)
out_DRT_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_eco_lai.SetGeoTransform(geotransform)
out_eco_lai.SetProjection(projection)
out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
out_eco_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_wat_lai.SetGeoTransform(geotransform)
out_wat_lai.SetProjection(projection)
out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
out_wat_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Float32)
out_tim_lai.SetGeoTransform(geotransform)
out_tim_lai.SetProjection(projection)
out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
out_tim_lai=None
print(rt_hv)Das obige ist der detaillierte Inhalt vonWie Python das GDAL-Modul verwendet, um Rasterdaten zu lesen und die angegebenen Daten zu filtern. Für weitere Informationen folgen Sie bitte anderen verwandten Artikeln auf der PHP chinesischen Website!
![[Web-Frontend] Node.js-Schnellstart](https://img.php.cn/upload/course/000/000/067/662b5d34ba7c0227.png)
