Python이 GDAL 모듈을 사용하여 래스터 데이터를 읽고 지정된 데이터를 필터링하는 방법
- 王林앞으로
- 2023-05-15 11:16:051213검색
1 코드 분할 설명
1.1 모듈 및 경로 준비
먼저, 사용되는 모듈과 래스터 이미지 저장을 위한 다양한 경로를 준비해야 합니다.
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_path는 자사 제품의 저장 경로, gl_file_path는 GLASS 제품의 저장 경로, out_file_path 두 래스터 간의 차이 처리 후 최종 결과의 저장 경로입니다. rt_file_path为自有产品的存放路径，gl_file_path为GLASS产品的存放路径，out_file_path为最终二者栅格做完差值处理后结果的存放路径。
1.2 栅格图像文件名读取与配对
接下来，需要将全部待处理的栅格图像用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.3 输出文件名称准备
前述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产品进行相减，因此配置了四个输出路径文件夹。
1.4 栅格文件数据与信息读取
接下来，利用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是将图层原有的缩放系数复原。
其次，上述代码第三段为获取栅格行、列数与投影变换信息。
1.5 差值计算与QA波段筛选
接下来，首先对自有产品与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()을 사용하여 처리할 모든 래스터 이미지를 가져오고 for를 사용해야 합니다. 루프 일괄 처리 작업을 준비하기 위한 루프입니다. 🎜 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)🎜그 중 이 글의 요건은 두 제품을 비교하는 것이기 때문에 먼저 두 제품의 hv 프레임 번호를 합쳐서 두 개의 원격탐사 이미지를 동일한 프레임 번호로 넣어야 합니다 따라서 자체 제품 파일 이름의 특성에 따라 .split()를 선택하여 문자열을 분할한 다음 원격 감지의 hv 프레임 번호를 가로챕니다. 영상. 🎜🎜1.3 출력 파일 이름 준비🎜🎜 출력 파일을 저장할 경로는 앞서 언급한 1.1 섹션에서 구성했지만 출력 파일 이름은 아직 구성되지 않았으므로 여기서 각각을 구성해야 합니다. error 최종 원격탐사 영상의 파일 저장 경로 및 이름입니다. 그 중 원격탐사 영상의 hv 번호를 출력 결과 파일명으로 직접 사용합니다. 🎜# -*- 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)🎜자사 제품의 LAI는 각각 4가지 알고리즘을 기반으로 구해지기 때문에 이 코드 부분은 4가지 부분으로 나누어져 있습니다. 차이점을 만들려면 각 알고리즘을 으로 비교해야 합니다. GLASS를 별도로 제품을 빼서 4개의 출력 경로 폴더가 구성됩니다. 🎜🎜1.4 래스터 파일 데이터 및 정보 읽기🎜🎜다음으로 gdal 모듈을 사용하여 두 종류의 래스터 파일인 .tif 및 .hdf를 읽습니다. > 읽을 그리드 이미지. 🎜rrreee🎜먼저 위 코드의 첫 번째 문단을 예로 들어보겠습니다. 그 중에서 gdal.Open()은 래스터 이미지를 읽고 .RasterCount는 래스터 이미지 밴드 수를 얻습니다. 래스터 이미지로 이미지의 각 밴드 정보를 배열 형식으로 읽어옵니다. 밴드 수가 1보다 크면 3차원이고 첫 번째 차원입니다. rt_raster_array[ 0]는 Array의 첫 번째 밴드를 취하는 것을 의미하며, 이 문서에서는 우리의 LAI 밴드입니다. product; rt_qa_array=rt_raster_array[1 ]는 두 번째 밴드를 꺼내는 것을 의미합니다. 이 문서에서는 우리 제품의 QA 밴드입니다. )는 래스터 이미지를 가져오는 것을 의미합니다. 첫 번째 밴드(여기서 일련 번호는 0가 아니라 1부터 시작함). where(rt_lai_array>30000,np.nan, rt_lai_array)는 np.where() 함수를 사용하여 첫 번째 밴드의 픽셀 >30000을 선택하는 것을 의미합니다. Array에서 nan으로 설정하고 다른 값은 변경하지 않습니다. 이 단계는 이미지에서 패딩 및 Nodata 값을 제거하는 것입니다. 마지막 문장 *0.001은 레이어의 원래 배율을 복원하는 것입니다. 🎜🎜두 번째로 위 코드의 세 번째 섹션은 래스터 행, 열 번호 및 투영 변환 정보를 얻는 것입니다. 🎜🎜1.5 차이 계산 및 QA 밴드 스크리닝🎜🎜다음으로 먼저 자사 제품과 GLASS 제품에 대해 차이 연산을 수행한 후 4개의 알고리즘을 별도로 추출해야 합니다. 🎜rrreee🎜그 중 위 코드의 처음 두 문장은 차이 계산과 데이터 반올림입니다. 데이터를 정수로 변환하면 결과 데이터 레이어의 데이터 양이 줄어듭니다(소수점을 저장할 필요가 없기 때문입니다). 🎜🎜그런 다음 QA 밴드를 기반으로 데이터 필터링 및 마스킹을 시작합니다. 실제로 다양한 원격 탐사 이미지(예: MODIS, Landsat 등)의 QA 밴드는 상대적으로 유사합니다. 이진 코드 문자열 이미지 품질, 정보 등으로, 서로 다른 비트가 종종 특성을 나타냅니다. 예를 들어, 아래 그림은 Landsat Collection 2 Level-2의 QA 밴드의 의미를 보여줍니다. 🎜🎜여기서 QA 밴드는 원래 십진수(일반적으로 원격탐사 이미지에서 공간을 절약하기 위해 QA 밴드는 십진수로 표기함)이므로 필요합니다. 그런 다음 필요한 바이너리 데이터 자릿수를 얻음으로써(이 문서에서는 당사 제품의 이 픽셀이 어떤 알고리즘에서 나오는지 결정할 수 있는 바이너리 자릿수입니다) , 이를 통해 이 픽셀의 값을 결정하여 각 알고리즘에 해당하는 픽셀을 함께 처리하도록 LAI 알고리즘을 구합니다. DRT_lai_dif_array를 포함한 4개의 변수는 각각 자체 알고리즘으로 얻은 픽셀을 제외한 4가지 알고리즘의 모든 픽셀을 나타냅니다. 이 방법을 선택한 이유는 나중에 마스크를 제거하면 바로 사용할 수 있기 때문입니다. 남은 것은 이 알고리즘 자체의 픽셀입니다. 🎜🎜그 중 위 코드에 주석 처리된 plt 관련 내용을 활용해 공간분포도를 그릴 수 있으니 관심 있으신 분들은 한번 사용해보시면 좋을 것 같습니다. 🎜1.6 结果栅格文件写入与保存
接下来，将我们完成上述差值计算与依据算法进行筛选后的图像保存。
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表示将其从内存空间中释放，完成写入与保存工作。
2 完整代码
本文所需完整代码如下：
# -*- 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)위 내용은 Python이 GDAL 모듈을 사용하여 래스터 데이터를 읽고 지정된 데이터를 필터링하는 방법의 상세 내용입니다. 자세한 내용은 PHP 중국어 웹사이트의 기타 관련 기사를 참조하세요!
성명：
이 기사는 yisu.com에서 복제됩니다. 침해가 있는 경우 admin@php.cn으로 문의하시기 바랍니다. 삭제