利用GDAL对中国范围长时间序列栅格数据的分析

1.数据为中国2000-2015年的ndvi数据，放置位置如下图。

Paste_Image.png

2.安装GDAL，主要通过下载.whl文件，通过pip安装。

3.对每个单元格的多年数据进行线性回归，并给出F值。

#coding=utf-8
import os
import os.path
import gdal
import sys
from gdalconst import *
from osgeo import gdal
import osr
import numpy as np
import mk
import pandas as pd
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

def Read(RasterFile):#读取每个图像的信息
    ds = gdal.Open(RasterFile,GA_ReadOnly)    
    if ds is None:
        print 'Cannot open ',RasterFile
        sys.exit(1)
    cols = ds.RasterXSize
    rows = ds.RasterYSize
    band = ds.GetRasterBand(1)
    data = band.ReadAsArray(0,0,cols,rows)  
    return data

def Readxy(RasterFile): #读取每个图像的坐标信息并返回     
    ds = gdal.Open(RasterFile,GA_ReadOnly)
    if ds is None:
        print 'Cannot open ',RasterFile
        sys.exit(1)
    cols = ds.RasterXSize
    rows = ds.RasterYSize
    band = ds.GetRasterBand(1)
    #data = band.ReadAsArray(0,0,cols,rows)
    noDataValue = band.GetNoDataValue()
    projection=ds.GetProjection()
    geotransform = ds.GetGeoTransform()
    return rows,cols,geotransform,projection,noDataValue

def WriteGTiffFile(filename, nRows, nCols, data,geotrans,proj, noDataValue, gdalType):#向磁盘写入结果文件
    format = "GTiff"   
    driver = gdal.GetDriverByName(format)
    ds = driver.Create(filename, nCols, nRows, 1, gdalType)
    ds.SetGeoTransform(geotrans)
    ds.SetProjection(proj)
    ds.GetRasterBand(1).SetNoDataValue(noDataValue)
    ds.GetRasterBand(1).WriteArray(data)    
    ds = None
rows,cols,geotransform,projection,noDataValue = Readxy("D:/china/2000.tif") 
print rows,cols
count=0
rootdir = "D:/china"
datalist=[]
for dirpath,filename,filenames in os.walk(rootdir):#遍历源文件
    for filename in filenames:
        if os.path.splitext(filename)[1] == '.tif':#判断是否为tif格式
             filepath = os.path.join(dirpath,filename)
             filedata1 = [[0.0]*cols]*rows
             
             filedata2 = np.array(filedata1)
             filedata3 = Read(filepath)
             count+=1
             datalist.append(filedata3)
             
dtarelist=[]
for m in range(len(datalist)):
    dtarelist.append((datalist[m].reshape((1,396, 697))))
#print dtarelist
#print dtarelist[2].shape
for a in range(len(datalist)):
  if a==0:
    dtaz=np.concatenate([dtarelist[a],dtarelist[a+1]],axis=0)
  if a>1:
    dtaz=np.concatenate([dtaz,dtarelist[a]],axis=0)
#print dtaz.shape

    
Arraytime=np.array([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])

Time=np.vstack([Arraytime,np.ones(len(Arraytime))]).T
temparray=[]
XL=np.zeros((396, 697),dtype=np.float)
XF=np.zeros((396, 697),dtype=np.float)
for r in range(396):
     for c in range(697):
         for k in range(16):
              temparray.append(dtaz[k,r,c])
         data={'Time':Time[:,0],'Arrayslope':temparray}
         df=pd.DataFrame(data)
         #Arrayslope=np.array(temparray)
         temparray=[]
         formula='Arrayslope~Time'
         model=ols(formula,df).fit()
         f=anova_lm(model)
         m=model.params.Time
         #m,b=np.linalg.lstsq(Time,Arrayslope)[0]
         XF[r,c]=f.iloc[0,3]
         XL[r,c]=m
     print r
        
WriteGTiffFile("D:/cal/NDVIMEAN_L.tif",rows,cols,XL,geotransform,projection,noDataValue, GDT_Float32)
WriteGTiffFile("D:/cal/NDVIMEAN_F.tif",rows,cols,XF,geotransform,projection,noDataValue, GDT_Float32)                

4.在上述代码的基础上，增加中值计算的函数和mk计算的函数。

···python
def median(x):
if len(x) < 1:
return None
else:
x.sort()
return x[len(x) // 2]
def mk_trend(x):
import math
import numpy as np
s=0
length=len(x)
for m in range(0,length-1):
for n in range(m+1,length):
if x[n]>x[m]:
s=s+1
elif x[n]==x[m]:
s=s+0
else:
s=s-1
#计算vars
vars=length(length-1)(2*length+5)/18
#计算zc
if s>0:
zc=(s-1)/math.sqrt(vars)
elif s==0:
zc=0
else:
zc=(s+1)/math.sqrt(vars)

#计算za    
zc1=abs(zc)
bl=0
if zc1>1.96:#对应95%的置信度
  bl=1
#计算倾斜度
    
ndash=length*(length-1)/2
slope1=np.zeros(ndash)
m=0
for k in range(0,length-1):
    for j  in range(k+1,length):
        slope1[m]=(x[j]-x[k])/(j-k)
        m=m+1
        
slope=median(slope1)
return slope,bl,zc

5.对网格循环部分修改代码，将mk计算嵌入网格循环中。
···python
temparray=[]
x_slope=np.zeros((396, 697),dtype=np.float)
x_bl=np.zeros((396, 697),dtype=np.float)
x_zc=np.zeros((396, 697),dtype=np.float)
for r in range(396):
     for c in range(697):
         for k in range(16):
              temparray.append(dtaz[k,r,c])
         #data={'Time':Time[:,0],'Arrayslope':temparray}
         #df=pd.DataFrame(data)
         #Arrayslope=np.array(temparray)
         
         #formula='Arrayslope~Time'
         #model=ols(formula,df).fit()
         slope,bl,zc=mk_trend(temparray)
         #f=anova_lm(model)
         #m=model.params.Time
         #m,b=np.linalg.lstsq(Time,Arrayslope)[0]
         x_slope[r,c]=slope
         x_bl[r,c]=bl
         x_zc[r,c]=zc
         temparray=[]
     print r
        
WriteGTiffFile("D:/cal/NDVIMEAN_slope2.tif",rows,cols,x_slope,geotransform,projection,noDataValue, GDT_Float32)
WriteGTiffFile("D:/cal/NDVIMEAN_bl2.tif",rows,cols,x_bl,geotransform,projection,noDataValue, GDT_Float32)
WriteGTiffFile("D:/cal/NDVIMEAN_zc2.tif",rows,cols,x_zc,geotransform,projection,noDataValue, GDT_Float32)  
···

其中x_slope为坡度值， x_zc为检验判断依据， x_bl为判断依据，值为1即通过mk检验。