CSV文件转矢量，并用克里金进行插值处理，矢量转栅格，python

# MakeXYLayer.py # Description: Creates an XY layer and exports it to a layer file import shapefile as shp import os import csv import glob # import system modules import arcpy from arcpy.sa import * # convert CSV to Lyr format files # Set environment settings #======================================================= arcpy.env.workspace = r"E:\Qixiang_data_processing\Each_hour_monitoring_data\Tem_CSV" input_dir = "E:\Qixiang_data_processing\Each_hour_monitoring_data" #====================================================== os.chdir(input_dir) #os.mkdir('Tem_SHP') #os.mkdir('Tem_TIF') #os.mkdir('Krigoutput') #======================================================= z_values = ['Max_tem','Min_tem','Mean_tem'] for value in z_values: #=================================================== file_list = glob.glob(input_dir + '\Tem_CSV\*.csv') for file in file_list: # Set the local variables in_Table = file x_coords = "lon" y_coords = "lat" z_coords = value basename = os.path.basename(file) name = basename.split('.',1)[0] out_Layer = name + '_' + value #=============================================== saved_Layer = os.path.join(input_dir+'\Tem_CSV',name)+'_'+value+'.lyr' # Set the spatial reference spRef = r"E:\Qixiang_data_processing\Each_hour_monitoring_data\prj\GCS_WGS_1984.prj" # Make the XY event layer... arcpy.MakeXYEventLayer_management(in_Table, x_coords, y_coords, out_Layer, spRef, z_coords) # Print the total rows # print(arcpy.GetCount_management(out_Layer)) # Save to a layer file arcpy.SaveToLayerFile_management(out_Layer, saved_Layer) # convert lyr files to shp format files #======================================================= lyr_list = glob.glob(input_dir+'\Tem_CSV\*.lyr') for lyr in lyr_list: inFeatures = lyr #=================================================== outLocation = input_dir + '\Tem_SHP' # Execute FeatureClassToGeodatabase arcpy.FeatureClassToShapefile_conversion(inFeatures, outLocation) # operate Kriging processing # Description: Interpolates a surface from points using kriging. # Set environment settings # Set local variables # the length of field must less than 10 characters # MakeXYLayer.py # Description: Creates an XY layer and exports it to a layer file #====================================================== arcpy.env.workspace = r"E:\Qixiang_data_processing\Each_hour_monitoring_data\Tem_TIF" search_dir = input_dir + '\Tem_SHP' shp_files = glob.glob(search_dir+'\*.shp') try: for inFeatures in shp_files: # the unit of cellsize is degree, so 1km is equal to 1/111, that is 0.009 print(inFeatures) cellSize = 0.01 #--------------------------------------------- raster_basename = os.path.basename(inFeatures) file_name = raster_basename.split('.',1)[0] name = (file_name+'.tif').replace(' ','') c = raster_basename.split('.',1)[0] d = c.index('_') field = c[d+1:][0:10] #--------------------------------------------- outRaster = name outVarRaster = file_name+'var.tif' kModel = "CIRCULAR" kRadius = 12 # Check out the ArcGIS Spatial Analyst extension license arcpy.CheckOutExtension("3D") # Execute Kriging arcpy.Kriging_3d(inFeatures, field, outRaster, kModel, cellSize, kRadius,outVarRaster) # Save the output # outKriging.save(input_dir+'\Krigoutput') except Exception as err: print(err)