目錄

• 1 簡介與技術流程
• 2 數據及其預處理
• 3 代碼
• • 3.1 要調用的包
  • 3.2 代碼函數：計算斜率和變異系數
  • 3.3 代碼：計算柵格圖像的斜率
  • 3.4 代碼：計算柵格圖像的變異系數
  • 3.5 代碼：函數的調用與結果
• 4 后記

1 簡介與技術流程

之前看到一篇文章，用變異系數（CV）計算年際初級凈生產力（NPP）的變異系數，研究一個地區的產量穩定性。下面是這篇文章（Hong C, Jin X, Ren J, et al. 2019. Satellite data indicates multidimensional variation of agricultural production in land consolidation area. Science of the Environment, 653: 735-747.）關于這部分的截圖，供大家參考。上面一張圖是他材料與方法中的介紹，下面一張圖是他的結果。


總的來說，就是把多幅影像疊加在一起，逐像元構建一組時間序列，然后計算它的變異系數。比方說，你研究區覆蓋的像元是100*100的，時間序列是從2011到2015每年一幅影像，那么你就需要構建10000組時間序列，每組里面包含5個數據。我的思想是，把它放在一個坐標系中，橫軸是時間，縱軸是每個時間下的像元值，一組時間序列構建一條折線，比較直觀的可以參考我論文里的圖，如下（我這篇論文還在撰寫，所以不能完全放出來，作了一些修改）。

基于這個思想，我們可以計算得出遙感柵格圖像的變異系數。然后我再發散一下，發現還可以通過這個計算變化度（斜率）。我結合的實例也是初級凈生產力（NPP），用斜率來探究一個地區的生產力增長情況（Production Improvement）；用變異系數來探究生產穩定性（Production Stability）。

2 數據及其預處理

首先，我通過CASA模型計算出4個年份下的初級凈生產力（NPP），柵格數據存成geotif格式（后綴是.tif）。

然后把tif文件放在一個找得到的位置，等下要用。

3 代碼

3.1 要調用的包

from osgeo import gdal import numpy as np from sklearn import linear_model import copy

osgeo里的gdal是用來處理tif數據的。
numpy是用來做數據分析統計什么的。
sklearn里的linear_model是用來用直線擬合時間序列求斜率的。
copy用里面的深拷貝以防破壞影像數據的原矩陣。

3.2 代碼函數：計算斜率和變異系數

首先，我們先寫一個普通的斜率和變異系數的函數做個鋪墊。這里寫的函數先不考慮處理影像，單獨就處理一個數組的斜率和變異系數。

# 計算斜率 def calculate_slope(data):reg = linear_model.LinearRegression()reg.fit(np.array(range(len(data))).reshape(-1, 1), np.array(data).reshape(-1, 1))slope = reg.coef_ # 斜率intercept = reg.intercept_ # 截距slope = slope[0][0]intercept = intercept[0]return slope# 計算變異系數 def coefficient_of_variation(data): # 變異系數mean = np.mean(data) # 計算平均值std = np.std(data, ddof=0) # 計算標準差cv = std/meanreturn cv

3.3 代碼：計算柵格圖像的斜率

然后開始涉及到柵格圖像的斜率計算。這里需要調用上面寫的函數。
這里輸入的參數images是tif數據路徑的數組形式，通過for in逐個讀取（要保證輸入的的多幅影像的空間范圍、像元數和坐標系一致）。outpath是計算后輸出的斜率圖像。

# 柵格圖像組計算斜率 def slope(images, outpath):images_pixels = [] # 存放多個圖像像元矩陣的空數組for image in images:print(f'正在讀取柵格圖像： {image}')tif = str(image)open_tif = gdal.Open(tif) # 打開柵格圖像band = open_tif.GetRasterBand(1).ReadAsArray() # 獲取波段的矩陣images_pixels.append(band) # 把圖像的矩陣加入到數組中SLP = copy.deepcopy(images_pixels[0]) # 獲取一個矩陣作為要寫入的模板# 讀取像元，并計算變異系數print('正在計算')for i in range(len(SLP)):for j in range(len(SLP[1])):if SLP[i][j] >= 0:SLP_data = [] # 存放ij坐標下像元值的數組，以計算變異系數for px in range(len(images)): # 遍歷多個圖像下ij坐標的像元值SLP_data.append(images_pixels[px][i][j]) # 同一坐標的多點加入數組，以計算變異系數SLP_value = calculate_slope(SLP_data) # 變異系數計算SLP[i][j] = SLP_value # 寫入該坐標下的變異系數print('計算完成')# 保存柵格圖像gtiff_driver = gdal.GetDriverByName('GTiff')out_tif = gtiff_driver.Create(outpath, SLP.shape[1], SLP.shape[0], 1, gdal.GDT_Float32)# 將數據坐標投影設置為原始坐標投影out_tif.SetProjection(open_tif.GetProjection())out_tif.SetGeoTransform(open_tif.GetGeoTransform())out_band = out_tif.GetRasterBand(1)out_band.WriteArray(SLP)out_band.FlushCache()print('柵格圖像組斜率計算完成')

3.4 代碼：計算柵格圖像的變異系數

與計算斜率類似。

# 柵格圖像組計算變異系數 def CV(images, outpath):images_pixels = [] # 存放多個圖像像元矩陣的空數組for image in images:tif = str(image)open_tif = gdal.Open(tif) # 打開柵格圖像band = open_tif.GetRasterBand(1).ReadAsArray() # 獲取波段的矩陣images_pixels.append(band) # 把圖像的矩陣加入到數組中CV = copy.deepcopy(images_pixels[-1]) # 獲取一個矩陣作為要寫入的模板# 讀取像元，并計算變異系數for i in range(len(CV)):for j in range(len(CV[1])):CV_data = [] # 存放ij坐標下像元值的數組，以計算變異系數for px in range(len(images)): # 遍歷多個圖像下ij坐標的像元值CV_data.append(images_pixels[px][i][j]) # 同一坐標的多點加入數組，以計算變異系數CV_value = coefficient_of_variation(CV_data) # 變異系數計算CV[i][j] = CV_value # 寫入該坐標下的變異系數# 保存柵格圖像gtiff_driver = gdal.GetDriverByName('GTiff')out_tif = gtiff_driver.Create(outpath, CV.shape[1], CV.shape[0], 1, gdal.GDT_Float32)# 將數據坐標投影設置為原始坐標投影out_tif.SetProjection(open_tif.GetProjection())out_tif.SetGeoTransform(open_tif.GetGeoTransform())out_band = out_tif.GetRasterBand(1)out_band.WriteArray(CV)out_band.FlushCache()print('柵格圖像組變異系數計算完成')

3.5 代碼：函數的調用與結果

我這里用的是2017年到2020年四年的四幅影像。影像長什么樣在“數據及其預處理”部分已經展示過了。
下面是執行處理的代碼。

images = [r'E:\2017.tif',r'E:\2018.tif',r'E:\2019.tif',r'E:\2020.tif', ]slope(images, r'E:\斜率2017_2020.tif') CV(images, r'E:\變異系數2017_2020.tif')

這里我把images單獨拿出來存成數組作為輸入，來保證代碼美觀（記得一定要按時間順序排）。
然后后面直接調用上面寫的slope函數和CV函數計算斜率和變異系數。images是輸入柵格時間序列圖像路徑，后面那個參數是輸出存放的路徑。
結果如下面，這里只展示斜率的。

我這里僅結合我自己的實驗做個展示，一個時間序列只有4幅影像，如果要做更長的幾十幅幾百幅的都不是問題，如果太多的話可以自己另寫個批量處理什么的。

4 后記

之前還是人文地理學專業的我，這兩個月才在學習遙感，可能有些地方不太專業不太科學，希望諸位同行前輩不吝賜教或者可以和我共同探討一下。可以在csdn私信我或者聯系我郵箱（chinshuuichi@qq.com），不過還是希望大家郵箱聯系我，我最近不經常上csdn。

如果覺得有幫助，希望大家能打賞支持一下（我的乞討碼在下面）~
-----------------------分割線（以下是乞討內容）-----------------------

總結

以上是生活随笔為你收集整理的【Python】基于Python计算长时间遥感栅格图像的像元值变化度（斜率）和变异系数的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。