Pertama sekali, anda perlu menyediakan modul yang digunakan dan pelbagai laluan untuk menyimpan imej raster.
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/"
Antaranya, rt_file_path ialah laluan penyimpanan produknya sendiri, gl_file_path ialah laluan penyimpanan produk KACA dan out_file_path ialah pemprosesan perbezaan terakhir antara dua raster. Laluan penyimpanan hasil akhir.
Seterusnya, semua imej raster yang akan diproses perlu diperolehi dengan os.listdir() dan gelung for digunakan untuk pemprosesan kelompok gelung Persediaan untuk operasi.
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 Antaranya, memandangkan keperluan artikel ini adalah untuk membandingkan dua produk, pertama sekali perlu menggabungkan hv nombor bingkai kedua-duanya dan meletakkan dua imej penderiaan jauh dengan nombor bingkai yang sama bersama-sama ; oleh itu, Berdasarkan ciri-ciri nama fail produk sendiri, pilih .split() untuk pembahagian rentetan, dan kemudian memintas nombor bingkai hv untuk mendapatkan imej penderiaan jauh.
Laluan untuk menyimpan fail output telah dikonfigurasikan dalam bahagian 1.1 yang disebutkan di atas, tetapi nama fail output belum lagi dikonfigurasikan; jadi di sini kita Adalah perlu untuk mengkonfigurasi laluan penyimpanan fail dan nama setiap imej penderiaan jauh selepas pemprosesan. Antaranya, kami terus menggunakan nombor hv imej penderiaan jauh sebagai nama fail hasil output.
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")Bahagian kod ini terbahagi kepada empat bahagian kerana LAI produk kami sendiri diperoleh berdasarkan empat algoritma masing-masing Apabila membuat perbezaan, setiap algoritma perlu digabungkan Produk GLASS ditolak, jadi empat folder laluan output dikonfigurasikan.
Seterusnya, gunakan modul gdal untuk membaca dua imej raster, .tif dan .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()Pertama, mari kita ambil perenggan pertama kod di atas sebagai contoh. Antaranya, gdal.Open() membaca imej raster; .RasterCount memperoleh bilangan jalur imej raster; lebih besar daripada .ReadAsArray(), yang Terdapat tiga dimensi, dan dimensi pertama ialah bilangan jalur Array bermakna mengambil jalur pertama dalam 1, yang dalam artikel ini ialah jalur rt_raster_array[0]LAIArray; produk kami sendiri; bermaksud Keluarkan jalur kedua, yang dalam artikel ini ialah jalur QArt_qa_array=rt_raster_array[1] produk kami sendiri bermaksud mendapatkan jalur pertama dalam imej raster (perhatikan bahawa siri nombor di sini tidak bermula dari bermula dari .GetRasterBand(1)); 0 bermaksud menggunakan fungsi 1 untuk memilih piksel np.where(rt_lai_array>30000,np.nan,rt_lai_array) dalam jalur pertama dalam np.where() dan tetapkannya kepada Array, meninggalkannya nilai lain tidak berubah. Langkah ini adalah untuk menghapuskan padding dan nilai >30000 dalam imej. Ayat terakhir nan adalah untuk memulihkan faktor penskalaan asal lapisan. Nodata*0.001Kedua, bahagian ketiga kod di atas adalah untuk mendapatkan baris raster, nombor lajur dan maklumat transformasi unjuran.
1.5 Pengiraan perbezaan dan saringan jalur QA
, kemudian empat algoritma perlu dilakukan secara berasingan ekstrak. 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()
Kemudian, saringan dan penyamaran data bermula berdasarkan jalur
QA. Malah, jalur QA pelbagai imej penderiaan jauh (seperti MODIS, Landsat, dsb.) adalah agak serupa: kualiti imej penderiaan jauh adalah diwakili oleh rentetan kod binari , maklumat, dsb., di mana bit yang berbeza sering mewakili ciri. Contohnya, rajah di bawah menunjukkan maksud kumpulan QA Koleksi Landsat 2 Tahap-2. Di sini, jalur
QA asalnya dalam perpuluhan (untuk menjimatkan ruang dalam imej penderiaan jauh umum, jalur QA ditulis dalam bentuk perpuluhan), jadi ia perlu ditukar kepada Perduaan; kemudian tentukan piksel dengan mendapatkan bilangan digit data perduaan yang diperlukan (dalam artikel ini, bilangan digit perduaan yang boleh menentukan algoritma mana piksel ini dalam produk kami sendiri berasal dari). Algoritma manakah yang digunakan untuk mendapatkan LAI yang diperoleh, supaya piksel yang sepadan dengan setiap algoritma diproses bersama. Empat pembolehubah seperti masing-masing mewakili semua piksel dalam empat algoritma kecuali piksel yang diperolehi oleh algoritma sendiri sebab memilih kaedah ini adalah kerana kita boleh menutupnya secara langsung kemudian, maka Yang tinggal hanyalah piksel; algoritma itu sendiri. DRT_lai_dif_array Antaranya, kandungan berkaitan
接下来，将我们完成上述差值计算与依据算法进行筛选后的图像保存。
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)Atas ialah kandungan terperinci Cara Python menggunakan modul GDAL untuk membaca data raster dan menapis data yang ditentukan. Untuk maklumat lanjut, sila ikut artikel berkaitan lain di laman web China PHP!
![[Web front-end] Permulaan pantas Node.js](https://img.php.cn/upload/course/000/000/067/662b5d34ba7c0227.png)
