安装

1、安装Anaconda，下一步

2、安装GDAL

（1）打开Anaconda prompt，输入conda install gdal
（2）打开网址：https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal
下载对应python的gdal包，注意对应的Python版本，如GDAL-3.2.3-cp37-cp37m-win_amd64对应的是Python37的，

（3）进入下载目录cd xxx
（4）执行pip install GDAL-3.2.3-cp37-cp37m-win_amd64.whl
（5）设置Windows环境变量
path:C:\Users\xxx\Anaconda3\Lib\site-packages\osgeo
GDAL_DATA:C:\Users\xxx\Anaconda3\Lib\site-packages\osgeo\data\gdal

（5）程序中使用：from osgeo import gdal

使用

简单的使用：打开影像，读取像素值，读取影像信息，保存影像等。

from datetime import datetime
from osgeo import gdal
import os  
import pandas as pd
import numpy as np
import date_util as du

threshold = 0.20
    
#  数组保存为tif
def array2raster(TifName, projection, geotransform, array):
    if os.path.exists(TifName):
        os.remove(TifName)
    cols = array.shape[1]  # 矩阵列数
    rows = array.shape[0]  # 矩阵行数
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(TifName, cols, rows, 1, gdal.GDT_Byte)
    # 括号中两个0表示起始像元的行列号从(0,0)开始
    outRaster.SetGeoTransform(geotransform)
    # 获取数据集第一个波段，是从1开始，不是从0开始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    #outRasterSRS = osr.SpatialReference()
    # 代码4326表示WGS84坐标
    #outRasterSRS.ImportFromEPSG(4326)
    #outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outRaster.SetProjection(projection)
    outband.FlushCache()

def ndvi_process(file_dir, output):
    file_names = os.listdir(file_dir)
    
    wid = 0
    hei = 0
    k = 0
    
    for file in file_names:
        ff = file_dir+"\\"+file
        print(ff)
        Raster1 = gdal.Open(ff)
        if wid == 0:
            wid = Raster1.RasterXSize
            hei = Raster1.RasterYSize
            
            #global resultarr, geotransform, projection
            resultarr = np.empty((hei, wid))
            
            geotransform = Raster1.GetGeoTransform()
            projection = Raster1.GetProjection()
            
        else:
            if wid != Raster1.RasterXSize:
                print("图像错误，大小不相等！")
                return
            
        data1 = Raster1.ReadAsArray()
        
        str = file[-12:-4]
        date_object = datetime.strptime(str, '%Y%m%d').date()
        val = ndvi_value(date_object)
        print(val)
        
        if val > 0:
            k = k + 1
            for i in range(hei):
                for j in range(wid):
                    # ndvi process
                    if data1[i][j] > (val-threshold) and data1[i][j] < (val+threshold):
                        resultarr[i][j] = resultarr[i][j] + 1
                    else:
                        resultarr[i][j] = 0

                
    array2raster(output, projection, geotransform, resultarr)

ndvi_process(r"D:\ndvi\ndvi_2018", r"D:\ndvi\area.tif")