Rumah > Artikel > pembangunan bahagian belakang > Cara Python menggunakan modul GDAL untuk membaca data raster dan menapis data yang ditentukan
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.001
Kedua, 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()
Antaranya, dua ayat pertama kod di atas ialah pengiraan perbezaan dan pembundaran data. Menukar data kepada integer mengurangkan jumlah data dalam lapisan data yang terhasil (kerana tidak perlu menyimpan perpuluhan).
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!