Arcpy: Spatial Join

Introduction

Points落在哪些Polygon之上，如Fig. 1？SpatialJoin_analysis解决之。

Fig. 1

注意，Points和Polygon文件属性表尽量简洁，留下关注的属性，不必要的属性在处理前删除。如Fig. 2和Fig. 3。

Fig. 2

Fig. 3

运行结束后，我们就可以了解到这些Points具体归属情况，通过属性表，见Fig. 4。

Fig. 4

Files.

References

[1] How to aggregate points to polygons using ArcGIS Desktop?

[2] Spatial Join (Analysis).

Arcpy: Merge

Merge Computes a geometric union of the input features. All features and their attributes will be written to the output feature class.

ArcToolbox--Data Management --General--Merge

注意文件使用绝对路径，否则可能引起不可预知的问题。Files

Reference

[1] Merge (Data Management), ArcGIS Help.

Python+Arcpy: Interpolating Using IDW or Kriging

已知：一个Point Shapefile文件，基于其中的字段完成由点及面的插值工作。

先以IDW反距离加权方法插值，方法：体会Idw，或者IDW_ga。插值如下：

Fig. 1

File.rar

Kriging未来再加入。

References

[1] How IDW works?

Create a Point Shapefile based on Lon/Lat imported by a Excel file

已知：Excel文件sheet表中包含气象站点的某类数据，并附上站点的经度与纬度，试探讨如何将该数据生成Point Shapefile文件，并为各个Point添加属性字段及赋值？

方法：Python+xlrd读取Excel文件中的sheet内容，Python+Arcpy实现Point Shapefile文件创建及其字段添加及赋值环节。

认真体会CopyFeatures_management、AddField_management、UpdateCursor。

运行结果如Fig. 1，红框表示添加字段及赋值后的属性表。File.rar.

Fig. 1

Matlab+Arcpy: 批量输出ArcMap图片

Introduction

配置要求tif与mxd同名，代码并行输出PNG图片，Layout View。

Fig. 1

File.rar

MXD匹配M代码.

禁止ArcMap自动拉伸 Go to Customize > ArcMap Options > Raster tab > Raster Layer sub tab and set the default to None.

References

[1] ExportToPNG, ArcGIS.

[2] Turn off automatic stretch when opening raster?.

Python+Arcpy: Points To Line

Introduction

Python+Arcpy操作Points（.shp）转换至Polyline（.shp），仔细研读Points To Line (Data Management)说明，参数说明如下：

Input_Features: The point features to be converted into lines.

Output_Feature_Class:The line feature class which will be created from the input points.

以下参数对生成闭合Polyline尤为重要（Optional）

Line_Field: Each feature in the output will be based on unique values in the Line Field.这里指出同一直线的各点应在某一Field下具有相同的数值，见Fig 1红框之内在Line1字段之下数值均为1，表示各点均在同一条Polyline。

Fig. 1

Sort_Field: By default, points used to create each output line feature will be used in the order they are found. If a different order is desired, specify a Sort Field.

Close_Line: Specifies whether output line features should be closed. True or False.

举一个例子，从含有4个point的.shp转换为1条闭合的Polyline（.shp）文件，转换输出如Fig 2。

Fig. 2

File.rar

注意，输入文件需使用绝对路径，起初使用相对路径返回错误指出无此文件。

Arcpy: Anselin Local Moran’s I

Introduction

Global spatial analysis or global spatial autocorrelation analysis yields only one statistic to summarize the study area. In other words, global analysis assumes homogeneity. If that assumption does not hold, then having only one statistic does not make sense as the statistic should offer over space.

But if there is no global autocorrelation or no clustering, we can still find clusters at a local level using local spatial autocorrelation. The fact that Moran's I is a summation of individual crossproducts is exploited by the "local indicators of spatial association" (LISA) to evaluate the clustering in those individual units by calculating Local Moran's I for each spatial unit and evaluating the statistical significance for each I_i, the Local Moran's I statistic of spatial association is given as:
\[{I_i} = \frac{{{x_i} - \overline x }}{{S_i^2}}\sum\limits_{j = 1,j \ne i}^n {{w_{i,j}}\left( {{x_j} - \overline x } \right)} \]

where x_i is an value of location i, xbar is the mean of the corresponding x, w_ij is the spatial weight between i and j, and

\[S_i^2 = \frac{{\sum\limits_{j = 1,j \ne i}^n {{{\left( {{x_j} - \overline x } \right)}^2}} }}{{n - 1}} - {\overline x ^2}\]

with n equating to the total number of pixels.

The z_Ii-score for the statistics are computed as:

\[{z_{{I_i}}} = \frac{{{I_i} - {\rm{E}}\left[ {{I_i}} \right]}}{{\sqrt {{\rm{V}}\left[ {{I_i}} \right]} }}\]

Where

\[{\rm{E}}\left[ {{I_i}} \right] =  - \frac{{\sum\limits_{j = 1,j \ne i}^n {{w_{ij}}} }}{{n - 1}}\]

\[{\rm{V}}\left[ {{I_i}} \right] = {\rm{E}}\left[ {I_i^2} \right] - E{\left[ {{I_i}} \right]^2}\]

\[{\rm{E}}\left[ {{I^2}} \right] = A - B\]

\[A = \frac{{\left( {n - {b_{{2_i}}}} \right)\sum\limits_{j = 1,j \ne i}^n {w_{i,j}^2} }}{{n - 1}}\]

\[B = \frac{{\left( {2{b_{{2_i}}} - n} \right)\sum\limits_{k = 1,k \ne i}^n {\sum\limits_{h = 1,h \ne i}^n {{w_{i,k}}{w_{i,h}}} } }}{{\left( {n - 1} \right)\left( {n - 2} \right)}}\]

\[{b_{{2_i}}} = \frac{{\sum\limits_{i = 1,i \ne j}^n {{{\left( {{x_i} - \overline x } \right)}^4}} }}{{{{\left( {{{\sum\limits_{i = 1,i \ne j}^n {\left( {{x_i} - \overline x } \right)} }^2}} \right)}^2}}}\]

The p-value for the null hypothesis is given as:

\[p = {\rm{erfc}}\left( {\frac{{\left| {{I_i} - {\rm{E}}\left[ {{I_i}} \right]} \right|}}{{\sqrt {2{\rm{V}}\left[ {I{}_i} \right]} }}} \right)\]

无奈，上面的这些公式计算结果就是不能与ArcGIS计算结果一致，只好放弃Matlab对栅格文件的直接操作，转而调用ArcPy操作。

Example

Fig.1

This tool creates a new Output Feature Class with the following attributes for each feature in the Input Feature Class: Local Moran's I index (LMiIndex), z-score (LMiZScore), p-value (LMiPValue), and cluster/outlier type (COType). The field names of these attributes are also derived tool output values for potential use in custom models and scripts. (GRID_CODE: Original Input Value)

Fig. 2

The z-scores and p-values are measures of statistical significance which tell you whether or not to reject the null hypothesis, feature by feature. In effect, they indicate whether the apparent similarity (a spatial clustering of either high or low values) or dissimilarity (a spatial outlier) is more pronounced than one would expect in a random distribution.

A high positive z-score for a feature indicates that the surrounding features have similar values (either high values or low values). The COType field in the Output Feature Class will be HH for a statistically significant (0.05 level) cluster of high values and LL for a statistically significant (0.05 level) cluster of low values.

A low negative z-score (for example, < -1.96) for a feature indicates a statistically significant (0.05 level) spatial outlier. The COType field in the Output Feature Class will indicate if the feature has a high value and is surrounded by features with low values (HL) or if the feature has a low value and is surrounded by features with high values (LH).

代码运行结果如Fig. 3所示，红框之处的COType显示非空，表示该点通过显著性检验（α=0.05）。

Fig. 3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

# Created by LI Xu # Version 1.0 # November 10, 2015 # Description: # Arcpy implemented Local Moran's I # If you have any question about this code, # please do not hesitate to contact me via E-mail: # jeremy456@163.com # Blog: # http://blog.sciencenet.cn/u/lixujeremy import arcpy import os import datetime import time def GetFileListExt(SouDir, FlagStr): import glob, os FileList=[] os.chdir(SouDir) filter_str='*'+FlagStr #print filter_str for filepath in glob.glob(filter_str): FileList.append(os.path.join(SouDir, filepath)) return FileList ##Main begintime=time.strftime("%Y-%m-%d %H:%M:%S") print 'Time Begin:'+begintime starttime = datetime.datetime.now() #Source Directory SouDir=r'E:\Tools\localMoranI\features' #Destination Directory DesDir=r'E:\Tools\localMoranI\features\output' #InputField InField='GRID_CODE' #Spatial Relationships SpRelations='INVERSE_DISTANCE' #Distance_Method DisMethod='EUCLIDEAN_DISTANCE' #Retrieve files FlagStr='.shp' files=GetFileListExt(SouDir,FlagStr) for filepath in files: filename=filepath.split('\\') filename=filename[len(filename)-1] filename=filename.split('.') otpath=DesDir+'\\'+filename[0]+'_LMI.shp' # Set geoprocessor object property to overwrite outputs if they already exist arcpy.gp.OverwriteOutput=True #Process: Local Moran's I clusters=arcpy.ClustersOutliers_stats(filepath, InField, otpath, SpRelations, DisMethod, 'NONE', '#', '#') print otpath print "END" endtime=time.strftime("%Y-%m-%d %H:%M:%S") print 'Time End:'+endtime finishtime=datetime.datetime.now() timespan=(finishtime-starttime).seconds timespan='%f' %timespan print 'Time Span:'+timespan+' s'