VIF代码实现

1.生成渔网

首先要把所有的因子统一分辨率。然后，按照栅格的分辨率来创建渔网。

下面代码是用python来统一分辨率

import arcpy
import os

# 设置工作空间
input_folder = r"E:\Patch-generating_Land_Use_Simulation_Model v1.4\workout_data\lulc_diver"
output_folder = r"E:\Patch-generating_Land_Use_Simulation_Model v1.4\workout_data\lulc_diver_re0.5k"
arcpy.env.workspace = input_folder

# 确保输出文件夹存在，如果不存在则创建
if not os.path.exists(output_folder):
    os.makedirs(output_folder)

# 列出工作空间中的所有栅格数据集
raster_list = arcpy.ListRasters()

# 循环处理每个栅格数据集
for raster in raster_list:
    # 获取栅格数据集的路径
    input_raster_path = os.path.join(arcpy.env.workspace, raster)
    
    # 设置输出分辨率为1000
    resolution = "500"
    
    # 使用 Resample_management 函数将栅格数据集转换为指定分辨率
    output_raster = os.path.join(output_folder, "resampled_" + raster)
    arcpy.Resample_management(input_raster_path, output_raster, resolution)

然后用ArcGIS中的创建渔网来创建渔网。
Data Management Tools-Sampling-Creat Fishnet

2.将值附在表格中

import arcpy
from arcpy import env
from arcpy.sa import *
import os




def ExtractValues_To_Points(filepath,inPointFeatures,inRaster,outPointFeatures ):
    """
    third factor
    NONE—不应用任何插值法；将使用像元中心值。 这是默认设置。
    INTERPOLATE—将使用双线性插值法根据相邻像元的有效值计算像元值。 将在插值中忽略 NoData 值，除非所有相邻像元均为 NoData。
    fourth factor
    VALUE_ONLY—仅将输入栅格的值添加到点属性。 这是默认设置。
    ALL—输入栅格的所有字段（Count 除外）都将添加到点属性。
    """

    # Set environment settings
    env.workspace = filepath
    arcpy.gp.overwriteOutput = True
    # Execute ExtractValuesToPoints
    ExtractValuesToPoints(inPointFeatures, inRaster, outPointFeatures,
                          "INTERPOLATE", "VALUE_ONLY")#INTERPOLATE


    # reference:https://pro.arcgis.com/zh-cn/pro-app/2.9/tool-reference/spatial-analyst/extract-values-to-points.htm


# v3版： 1、需要改进算法，不要重复输入


i = 1
while i < 12:

    Rasterpath = r"E:\Patch-generating_Land_Use_Simulation_Model v1.4\workout_data\lulc_diver_re1k"

    ShpPath =r"ES_1km_point{}.shp".format(i) #通过不断生成新的值来达到目的

    rasterList = arcpy.ListRasters("*", "TIF")

    ExtractValues_To_Points(Rasterpath,ShpPath,rasterList[i-1].format(i+1),"ES_1km_point{}.shp".format(i+1))

    # Set local variables
    inFeatures = "ES_1km_point{}.shp".format(i+1)
    fieldName = "bio{}".format(i + 1)
    arcpy.management.AddField(inFeatures, fieldName, "FLOAT")


    expression = "!RASTERVALU!"

    arcpy.management.CalculateField(inFeatures, fieldName, expression, "PYTHON3")

    # 参考：file:///C:/Program%20Files/ArcGIS/Pro/Resources/help/pro-app/3.0/tool-reference/data-management/calculate-field.htm


    # delect field

    fields = ["RASTERVALU"]
    method = "DELETE_FIELDS" #删除指定的值


    # Execute DeleteField
    arcpy.DeleteField_management(inFeatures, fields, method)
    print("完成——" + str(rasterList[i-1]))
    i= i+1

后面将矢量数据的栅格表导出成dbf，然后再导出CSV

3.R语言中做分析

在R语言中做分析

library(dplyr)
library(ggplot2)
library(caret)
library(car)
library(tidyverse)

getwd()

incsv <- read.csv("E:\\Patch-generating_Land_Use_Simulation_Model v1.4\\workout_data\\lulc_diver_re1k\\all_factor.csv")


indf <- data.frame(incsv)


summary(indf) # 得到特征值的

indf[indf=='-9999'] <- NA
indf=na.omit(indf)



head(indf)


indf_v1 <- subset(indf, select = -c(Id,bio2,bio8))


lm_v1 <-lm(indf_v1)

alias(lm_v1)

vif(lm_v1)