FVC值的计算

# coding=utf-8

import pandas as pd
import numpy as np
import os
import arcpy
from arcpy.sa import *


def load_img_to_array(img_file_path):
    inRas = arcpy.Raster(img_file_path)
    lowerLeft = arcpy.Point(inRas.extent.XMin, inRas.extent.YMin)
    ndarray = arcpy.RasterToNumPyArray(inRas)
    ndarray = ndarray.reshape(1, -1)
    ndarray.sort()
    ndarray = ndarray[ndarray>=-1] # 取正常的NDVI值
    q_5 = np.percentile(ndarray, 5)
    q_15 = np.percentile(ndarray, 15)
    q_30 = np.percentile(ndarray, 30)
    q_60 = np.percentile(ndarray, 60)

def get_confidence(filepath):
    inRas = arcpy.Raster(filepath)
    lowerLeft = arcpy.Point(inRas.extent.XMin, inRas.extent.YMin)
    ndarray = arcpy.RasterToNumPyArray(inRas)
    ndarray = ndarray.reshape(1, -1)
    ndarray.sort()
    ndarray = ndarray[ndarray >= -1]  # 取正常的NDVI值
    q_5 = np.percentile(ndarray, 5)
    q_95 = np.percentile(ndarray, 95)

    return q_5, q_95


arcpy.CheckOutExtension("spatial") # 检查外部扩展
arcpy.gp.overwriteOutput = 1 # 覆盖之前的文件

dirs = r'D:\dataset\caijiannc\ndvi/'
out_dir = r'D:\dataset\caijiannc\FVC/'
files = os.listdir(dirs)
files = [x for x in files if ('tif' in x)&(len(x) == 14)]
print(files)
for file in files:
    if os.path.exists(out_dir + file ) == False:
        os.mkdir(out_dir + file)
    out_path = out_dir + file + os.sep
    print(out_path)
    dir = dirs + file + os.sep
    print(dir)
    arcpy.env.workspace = dir
    rasters = arcpy.ListRasters(raster_type='TIF')

    for raster in rasters:
        inRaster = Raster(raster)
        Min, Max = get_confidence(dirs + file)
        print(Min, Max)
        ans = Con(inRaster < Min, 0, Con((inRaster >= Min) & (inRaster <= Max), (inRaster - Min)/(Max - Min), 1))
        ans.save(out_path + file)
        print(out_path + file)
        print (file + " is done!")

print ("OK!")