Zhuotong Nan (南卓铜, zhn1@pitt.edu)

DEM processing

In this part, we will show you how to delineate a watershed with inputs of DEM(Digital elevation model) and outlet location. Watershed delineation is necessary for any distributed hydrologic model. Although hydrologic models use different software to implement the delineation, the underlying idea and steps are generally common.

1. Open PIHMgis by double clicking the PIHMgis_v2.0.0 shortcut created in the previous section.

2. Enable PIHMgis.

Open Plugin Manager… in the QGIS Plugins menu, tick the PIHMgis item (Figure 3), and press Ok.

Figure 3 QGIS Plugin Manager

Now the PIHMgis toolbar appears as shown in Figure 4.

Figure 4 the PIHMgis v2.0.0 toolbar

3. Load DEM to PIHMgis for displaying.

Layer > Add a Raster Layer…, select lab3.shale.hills\DEM\sh_elev.asc.

This DEM file is in ESRI ASCII grid file format. If you have files in ESRI binary grid format which is a proprietary format, you need to convert them to ASCII format by using some ArcInfo command or by some third party tool, such as gdal_translate available from the GDAL website at http://www.gdal.org/.

The map view of the DEM file looks like as following (Figure 5),

Figure 5 DEM map view in QGIS

4. Save as a QGIS project. File > Save Project As…, in the Save As dialog, navigate to the directory of lab3.shale.hills, type “sh” for the file name, and then click Save.

Note: During following steps, please save the project frequently. If the project is crashed, you can restore to the previous session by open the project you saved.

5. Fill pits.

Question: Why do we need to fill pits for a DEM file prior to further watershed delineation?

Click the small arrow next to the Raster Processing on the PIHMgis toolbar and then select Fill Pits. This will bring up the Fill Pits dialog (Figure 6).

Figure 6 The Fill Pits interface

Select sh_elev.asc as the input DEM grid and name the output pit-filled grid as “pits_filled.asc” under the “raster_processing” directory as shown in Figure 6. Press Run to start the pit filling process. Click “Close” to dismiss the dialog after the process is done. The pit-filled DEM will automatically be loaded into the QGIS data frame.

6. Calculate flow direction and flow accumulation.

Use the command Raster_Processing > Flow Grid to open the Flow Grid window (Figure 7).

Set pits_filled.asc as the input; flow_dir.asc and flow_accu as the output flow direction grid and flow accumulation grid respectively, also under the raster_processing folder.

Click Run to start the processing. After done, click Close.

Figure 7 Flow Grid

Both flow_accu and flow_dir are shown in the QGIS map view (Figure 8).

Note: the direction codes are following: 1 - east, 2 - North east, 3 - North, 4 - North west, 5 -West, 6 - South west, 7 - South, 8 - South east.

Figure 8 Flow direction generated by PIHMgis

Tip: Within the QGIS framework, visibility of a GIS layer can be switched using the table of contents on left.

Tip: The identify tool can be used to identify features by clicking on them in the right map view panel. For example, click on the flow_dir item in the table of contents to make it selected, and then use the identify tool to click over the map view, the Identify Results dialog will pop up and show the value of the clicked feature (Figure 9), here, value of 4 means the point of interest has a flow direction towards north west.

Figure 9 Identify results

7. Generate stream grid

Select Raster Processing > Stream Grid to make the Stream Grid dialog visible, set the flow accumulation grid generated in previous step, i.e., flow_accu.asc, as the input grid, and then set output stream grid to stream_800.asc with a threshold value of 800 (number of grids). Press Run to generate stream grid (Figure 10).

Figure 10 Generate stream grid with a threshold value of 50.

Don’t close the Stream Grid window. Set thresholds to 300 and 800 respectively to compare the resulting stream grids (stream_300 and stream_800 accordingly, Figure 11) and run again.

Figure 11 The generated stream grids with thresholds set to 50, 300 and 800. (left: 50, middle: 300, right: 800)

Close the Stream Grid dialogue by clicking on the Close button.

Question: How to determine an appropriate threshold for a watershed of interest?

8. Generate the stream order grid (known as Link Grid in PIHMgis)

Select Raster Processing > Link Grid to pop up the configure interface (Figure 12). Select stream_800.asc as the input stream grid and flow_dir.asc as the input flow direction grid, and set the output link grid name to link_800.asc (saving to the same raster_processing folder). Click Run to begin processing. After done, the link grid will be shown on the map view panel if the Load in Data Frame option is checked.

Figure 12 Generate a stream order grid

9. Create a vector file of stream network

Select Raster Processing > Stream PolyLine to make the interface appear (Figure 13). Set inputs with stream_800.asc and flow_dir.asc generated by previous steps and set the name of the output stream polyline file to str800.shp. Click Run to process it after parameters are specified.

Figure 13 Create a vectorized stream network

10. Generate catchment grid and its corresponding vector file

Select Raster Processing > Catchment Grid to bring up the catchment grid generation interface (Figure 14). Set Stream Link Grid to link_800.asc, Flow Dir. Grid to flow_dir.asc, and output Catchment Grid to “cat_800.asc”. Both input and output files are located in the raster_processing directory.

Click Run to start calculating the catchment grid. Click Close to dismiss the dialog.

Select Raster Processing > Catchment Polygon to display the interface for converting grid catchment file to vector (Figure 15). Set input to cat_800.asc and output to cat800v.shp. Run to start the processing and Close to dismiss the window. The resulting catchment polygons are shown on Figure 16. Note: the color to render the polygons might not be same as that on the Figure 16, as it is determined by QGIS automatically.

Figure 14 Catchment grid generation interface

Figure 15 Convert catchment grid to vector

Figure 16 Resulting catchment polygons

11. Extract the watershed boundary of interest

The current PIHMgis version 2.0 does not provide such a function to clip the watershed boundary under study and stream network within it. We need to do this by help of other GIS software such as ESRI ArcInfo, or by using QGIS editing capability demonstrated below. In the future version of PIHMgis, this functionality will be implemented. To determine the boundary of a watershed, the outlet location should also be specified.

Open QGIS 1.0, load cat800v.shp and str800.shp generated in previous steps to the data frame using Layer > Add Vector Layer…, and then load outlet.shp, which contains outlet location. The view now looks like Figure 17.

Figure 17 QGIS 1.0 with stream network, catchment polygons and outlet layers loaded

Click on cat800v within the table of contents to make it the current layer. Use “Select Features” tool (Figure 18) to select the watershed polygon of interest (Figure 19). The selected polygon is now in yellow. If the tool is invisible, just drag the toolbar to a new place to make fully appeared.

Students are advised to have a look at all built-in toolbars, making them familiar with all those often-use tools.

Figure 18 The Attributes toolbar

Figure 19 Select the watershed of interest

Right click the cat800v item and select “Save selection as shapefile…” to save the selected polygon (Figure 20) to a new shapefile, here, “cat800v_sel.shp”, under the “raster_processing” folder. If you are required to select a project, click OK to use the default one.

Note: The default project probably is not correct, but does not matter with following operations. You can change the project in a layer’s property window later.

Figure 20 Save selection as a new shapefile

Similarly, select the corresponding stream from str800 and save as a new shapefile “str800_sel.shp” (Figure 21), also into the “raster_processing” folder. If you are required to select a project, click OK to use the default one.

Figure 21 Select the corresponding stream, shown in yellow.

Load str800_sel and cat800v_sel into the QGIS map view. Drag layers up or down to rearrange layers’ visibility. Uncheck str800 and cat800v to hide them. Now the map view looks like Figure 22.

Figure 22 Load str800_sel and cat800v_sel and hide str800 and cat800v

Make str800_sel disappear. Select “cat800v_sel” layer, and use “zoom in” tool (the one w/ magnifier icon and a plus symbol shown in Figure 23) to enlarge the area near the outlet. Click “Toggle editing” tool (the one w/ pen icon shown in Figure 23) to make the current layer editable (Figure 23). Now the other tools on the Digitizing toolbar are also available when the layer enters into edit mode. The tools may be hidden, in this case drag toolbars to a blank place to make all tools shown.

Figure 23 cat800v_sel in edit mode

Click the Split Features tool. Draw a curve to cross the outlet so as to split the whole polygon. Right click to finalize the curve; the Split Features will then split the polygon. Click the “Toggle editing” again to make the layer back in the normal view mode. Make sure saving it if you are asked to save the features.

Use “Select Features” tool to select the large polygon right to the outlet. If the previous operation has been done successfully, only the large part will be selected and the upper left part remains unselected. Save the selected polygon to a new shapefile named “cat800v_clip.shp”. Load cat800v_clip to the map view. It now looks like as shown on Figure 24.

Figure 24 Clipped watershed boundary

Similarly, we also split the str800_sel.shp. Make it visible and editable, draw a curve to split it across the outlet, save the operation, and then export the desirable line (on the right) to a new shapefile, “str800_clip.shp”. (Figure 25)

Figure 25 Clipped stream network

Zoom in the outlet area, and carefully edit the stream polyline making its end node locate on the watershed edge. However no necessary exact alignment is needed. Tools for vertex editing are shown in Figure 26. The final stream line looks like as in the Figure 27. Save the stream file after done.

Figure 26 vertex editing tools. Left: move vertex; middle: add vertex; right: remove vertex

Figure 27 Align the end node of stream line to the boundary. Left: original; right: edited

Close QGIS 1.0.

Zhuotong Nan (南卓铜, zhn1@pitt.edu)

Set up work place

1. Download PIHMgis version 2.0 from the URL:

http://www.pihm.psu.edu/PIHMgis_v2.0.0/PIHMgis_v2.0.0.zip

Download Shale Hills PIHMgis data files from the URL:

http://cid-0ea641a5a7f665a1.skydrive.live.com/self.aspx/Public/shale.hills.pihmgis.data.zip

Download qgis.1.0 from the URL:

http://cid-0ea641a5a7f665a1.skydrive.live.com/self.aspx/Public/qgis.1.0.zip

Note: The current PIHMgis v2.0 is built with QGIS 0.8 preview version and cannot run with QGIS 1.0 because APIs of those two version are inconsistent. The stable QGIS 1.0 version provides us a much better user interface and richer functionality than version 0.8. We will use QGIS 1.0 editing functions to manipulate vector data. The development team of PIHMgis is working on the migration to QGIS 1.0, which might be available in near future.

2. Make a directory “PIHM-lab”, and then extract “PIHMgis_v2.0.0.zip” to this directory. Go to this PIHMgis_v2.0.0 folder, locate “PIHMgis_v2.0.0.exe”, make a shortcut and rename the shortcut to “PIHMgis_v2.0.0”, and then copy this shortcut file to the directory “PIHM-lab” (our lab root directory).

Clicking the shortcut will bring up the Quantum GIS (QGIS) window, as shown Figure 1. Click “X” to close QGIS.

Figure 1 the QGIS 0.8 preview version interface coming with the PIHMgis. Note, the interface might be shown in different language if your computer locale is not set to English.

3. Create a directory “lab3.shale.hills” under “PIHM-lab”. Extract “shale.hills.pihmgis.data.zip” to the “lab3.shale.hills” directory.

4. Extract qgis.1.0.zip to the directory “PIHM-lab”, go to “qgis.1.0\bin”, make a shortcut “qgis” for qgis.exe under this folder, and then copy this shortcut to the directory “PIHM-lab”. Click the qgis shortcut to open it, making sure it can work. It looks like Figure 2.

Figure 2 the QGIS 1.0 interface

南卓铜 (zhn1@pitt.edu)

一些有名的Land cover数据集并没有使用通用的GIS格式，比如在以下链接，

http://www.geog.umd.edu/landcover/tree-cover-poster/download.html

下载broadleaf Layer，gzip解压缩后为 694,417,920 bytes，扩展名为.img，但不同于Erdas Imagine .img格式。

据介绍Goode binary file是plain binary，即相当于BSQ （Band sequential），其metadata有以下信息：

40031 pixels by 17347 lines，且是8bit unsigned integer，应当总大小为 694,417,757bytes，而解压缩后的大小大了163bytes，目前我还不知道原因。

根据metadata信息，可以添加头文件，broadleaf.hdr，内容如下，

nrows      17347
ncols     40031
nbands     1
nbits     8
layout     bsq
skipbytes     0
ulxmap     -20015500.000
ulymap     8673500.000
xdim     1000
ydim     1000

将broadleaf.img改名为 broadleaf.bsq，在ArcGIS中才可以正常访问，也可以导出为其它格式，比如Grid。不过注意其空间参考系仍为unknown。在ArcCatalog中打开其属性页，在spatial reference中选projected coordinates/World/Goode Homolosine (Land).prj。可以将bsq设置为 Goode投影。

如果希望转换成其它常用格式，

上图将bsq转成ESRI binary grid格式。图形化结果如下图。

在Workstation下好像没有找到Goode投影，但在Desktop下可以找到.prj。因为如果bsq是unknown时，导成esri grid，在ArcCatalog里对grid进行定义投影（属性->define spatial reference），没法子定义到Goode投影，很奇怪。所以这里，我是先对bsq定义投影，再转成grid。此外，通过ArcToolbox的Data Management下的projection define来做，也是可以。

When I run Identity_analysis in ArcGIS command windows against geodatabase, a message of something like "Number of shapes does not match the number of table records" popups. The error continues even if you reinstall the ArcGIS whenever using the Identify command.

Solution: go to c:\documents and settings\yourname\local settings\temp and delete all the files (some files might be denied to access, just skip them).

Zhuotong Nan (zhn1@pitt.edu)

以下简单的代码实现了ArcGIS支持的栅格向Golden software的Surfer ASCII grid的格式转化。需要注意的是，ArcGIS的栅格左上为坐标原点，向下向右为正。而Surfer grid是左下是坐标原点，向上向右为正。

Here I will show a class that implements the major functionality to convert any Raster file supported by ArcGIS to Surfer ASCII grid format. The resulting Surfer grid can be read by Surfer with version 6 or higher. One point we need to keep in mind, ArcGIS raster stores data in row-major order, right and downwards being positive, while Surfer grid takes lower left corner as the original, right and upwards being positive.

using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.Geometry;

namespace RasterToSurferGrid
{
    /// <summary>
    /// A wrapper class for Surfer ASCII grid file
    /// </summary>
    class SurferGrid
    {
        string id = "DSAA";
        int nx;//columns
        int ny; //rows
        double xlo; //the minimum X value of the grid
        double xhi; //xhi is the maximum X value of the grid
        double ylo;
        double yhi;
        double zlo;
        double zhi;
        List<double> data; //the grid is stored in row-major order, with the lowest row (minimum Y) first

        public SurferGrid(string rasterPath)
        {
            data = new List<double>(nx * ny);

            IWorkspaceFactory wf = new RasterWorkspaceFactoryClass();
            IRasterWorkspace wk=(IRasterWorkspace)wf.OpenFromFile(System.IO.Path.GetDirectoryName(rasterPath), 0);
            IRasterDataset ird = wk.OpenRasterDataset(System.IO.Path.GetFileName(rasterPath));
            IRaster raster=ird.CreateDefaultRaster();

            IRasterProps props = (IRasterProps)raster;
            nx = props.Width;
            ny = props.Height;
            IEnvelope env = props.Extent;
            xlo = env.XMin;
            xhi = env.XMax;
            ylo = env.YMin;
            yhi = env.YMax;
            //zlo = env.ZMin;
            //zhi = env.ZMax;

            //populate data
            PntClass extent=new PntClass();
            extent.X=nx;
            extent.Y=ny;
            IPixelBlock datablock=raster.CreatePixelBlock(extent);
            PntClass orig = new PntClass();
            orig.X = 0;
            orig.Y = 0;
            raster.Read(orig, datablock);

            for (int y = ny - 1; y >= 0; y--)
            //for (int x = 0; x < nx; x++)
            {
                for (int x = 0; x < nx; x++)
                {
                    data.Add((float)datablock.GetVal(0,x,y));
                }
            }

            //get the max, and min
            double max, min;
            max = data[0]; min = data[0];
            for (int i = 0; i < data.Count; i++)
            {
                max = max > data[i] ? max : data[i];
                min = min < data[i] ? min : data[i];
            }

            zlo = min;
            zhi = max;
        }

        public void WriteToFile(string outPath)
        {
            //StringWriter sw = new StringWriter();
            StreamWriter sw = new StreamWriter(outPath);
            sw.WriteLine(id);
            sw.WriteLine("{0} {1}", nx, ny);
            sw.WriteLine("{0} {1}", xlo, xhi);
            sw.WriteLine("{0} {1}", ylo, yhi);
            sw.WriteLine("{0} {1}", zlo, zhi);
            StringBuilder sb=new StringBuilder();
            for (int i = 0; i < ny; i++)
            {
                for (int j = 0; j < nx; j++)
                {
                    sb.AppendFormat("{0} ", data[i*nx+j]);
               }
               sb.Remove(sb.Length - 1, 1);
               sb.AppendLine();

            }
            sw.Write(sb.ToString());
            sw.Close();
        }

    }

}

附件是源代码。工具的运行需要dotnet framework 2.0和ArcView以上的License运行。编译后在命令行敲入RasterToSurferGrid有使用提示。

The attached please find necessary source codes which can be compiled on your own platform. Dotnet framework 2.0 and a valid ArcView license or higher are both required to run this tool. Besides, ArcGIS assemblies for dotnet framework are required. After successful compilation, typing RasterToSurferGrid following the command prompt shows you its usage information.

Zhuotong Nan(南卓铜, zhn1@pitt.edu)

教科书上都有MODIS有固定的过境时间，比如Terra说是10am过境（指从北向南经过赤道，中文资料上有意无意忽略这些定语，时间是UTC），Aqua是13:30过境（从南到北经过赤道。国家MODIS网站上将13:30写成14:30，是错误的。）。其实对于某一区域，过境时间不一样了。比如下面两张NASA Orbit Archive的Terra轨道图。

 
Aug 2, 2002(左），Dec 8, 2005(右)

对于美国Oklahoma区域(MODIS tile: h09v05)左边在下午5点，右边在下午7点。

而针对16天合成的NDVI和EVI，它是选择其后16天内的观测良好（如无云）的若干观测通过一定算法合成的，每个网格上的数据都有可能来自此16天内的不同时间。幸好NDVI/EVI作为反应植被的一种指标，对半月期变化是可以接受的，用来反应整个植被变化的趋势是足够了。但如果要跟土壤水分等时间性十分强（如一场大降水）的数据耦合到一起应用，则需要慎重的考虑。

(最近看的资料有这些想法，了解不多，如有谬误，请指正)

Zhuotong Nan (南卓铜，zhn1@pitt.edu)

网上可以找到的计算Kappa系数的工具都是基于ArcView 3.x，如KappaStat。大多是做Land cover classification的人开发的。计算对象是采样的点文件和分类的polygon或者grid文件。但如果我们将Kappa应用于两个图像比较它们的相似性时，这类软件不能直接用。一个可行的方法是，将其中一个图像转换成point文件，比如应用arcgis的RasterToPoint命令。但对于很多图像对要比较，则没有好用的法子。KappaStat的图形界面和本身是一个ArcView 3.x的扩展模块，决定了即使写Avenue脚本也不是件容易实现的事。

所以我们写了一个小工具Kappa.exe。它用来计算两个栅格文件之间的Kappa系数，只要两个栅格文件有同样的width和height（即两图像是同样的mxn阵列）。Kappa工具不考虑栅格文件所在的空间参考系，只是简单的逐网格进行运算。用户需要自己确认两幅图像有合适的空间参考系。

Kappa工具与KappaStat进行了结果验证，证明计算是可靠，输出结果包括误差矩阵（Error Matrix），Kappa系数，还包括Variance和P统计值。

Kappa工具对栅格图像的支持调用了ArcGIS ArcObjects栅格API，支持很多种栅格格式，包括ESRI grid, ASCII grid, TIFF/GeoTIFF, JPEG等。

Kappa工具运行在命令行。这意味着通过简单的批处理脚本，可以实现批处理功能。比如，以下Windows batch代码实现了两个目录下相同文件名的栅格文件的Kappa计算。计算结果保存在kappa.txt文件里。

@echo off
for %%i in (00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23) do (
echo Hour %%i ...
kappa.exe data1\rldas%%i data2\rldas%%i  -k:.\kappa.txt
)
echo Done.

命令行格式，Kappa.exe [参考或实测图像] [待比较图像] {-l:日志文件} {-k:Kappa文件}

[]表示必须，{}表示可选。对ESRI grid，图像指定到目录名；对个人geodatabase（通过Access数据库保存）的Raster dataset，以geodb.mdb\raster指定；文件geodatabase的Raster dataset，以geodb.gdb\raster指定。如果指定日志文件，在屏幕上输出的内容将同时保存到日志文件中；如果指定Kappa文件，则将Kappa系数单独保存在此文件中。此两选项开关在批处理有用。

运行需求

dotnet framework 2.0
ArcGIS desktop 9.2 with ArcView License or higher

运行平台

运行在Windows系统。但由于工具由C#编写，调用ArcGIS ArcObjects相关API，如果在Linux上有安装ArcGIS desktop 9.2，结合MONO runtime，本工具应当可以被编译成适合于Linux平台。

Binaries和source codes可以Email联系我。

Zhuotong Nan（南卓铜） zhn1@pitt.edu

ArcGIS 9.2提供了脚本功能，可以方便地实现批处理功能。可用脚本包括JS, VBS和Python。以下示例利用Python实现批处理ASCII grid文件到ESRI Binary grid转换的功能。

问题

ASCII grid由matlab代码生成，存储于 LDAS.data 目录。ASCII grid内容如下：

ncols         32
nrows         32
xllcorner     -88
yllcorner     37.75
cellsize      0.125
NODATA_value  -9999
0.000000    0.000000 ...

各ASCII grids命名有以下规律：

nbwi3.[yyyy].[mm].[dd].[hh]h.txt

其中[yyyy]是4位年份，[mm]是月份，[dd]是2位日期，[hh]是24制的小时。比如，nbwi3.2001.07.08.15h.txt等。

生成的ESRI Binary grid有以下命名规律：n[yy][mm][dd][hh]，其中[yy]是2位年份，[mm]是月份，[dd]是2位日期，[hh]是24制的小时。如n01070814。注意，ESRI grid的命名和长度是有一定要求限制的，所以我们必须要选择合适的，并能区别数据的命名方式。

我们需要做的是将数百个ASCII grids转换成对应的ESRI binary grids。

实现

脚本编写环境我们这里采用免费的PyDev/Eclipse。Python脚本存储在与LDAS.data并列的src目录下。事实上，可以使用任何文本编辑软件来用脚本，但Eclipse提供很多编程增强功能，如智能提示等。

import arcgisscripting
import os

datadir=os.path.abspath(os.path.join(os.getcwd(),"../LDAS.data"))
#print datadir
children=os.listdir(datadir)

#only ".txt" files
txtfiles=[]
grdnames=[]
for child in children:
    root,ext=os.path.splitext(child)
    if (ext=='.txt'):
        txtfiles.append(child);
        grdname='n'+root[8:19]
        grdname=grdname.replace('.','')
        grdnames.append(grdname)
# grid names now ready       

gp=arcgisscripting.create()
gp.workspace=datadir

# begin iteration
for i in range(len(grdnames)):
    gp.ASCIIToRaster(txtfiles[i], grdnames[i], "FLOAT")
    print txtfiles[i]

首先，需要通过import arcgisscripting，将arcgisscripting导入Python。arcgisscripting事实上是位于C:\Program Files\arcgis\bin目录下的arcgisscripting.dll。如果在运行时，提示找不到arcgisscripting库，则需要配置Python库路径。

库导入成功后，gp=arcgisscripting.create() 创建script对象，此后通过“gp.命令”的方式调用script支持的命令。一般来讲，可以在arcgis command window里使用的命令，都会有对象的script调用方式。比如，

ASCIIToRaster的命令语法是ASCIIToRaster_conversion <in_ascii_file> <out_raster> {INTEGER | FLOAT}，对应的脚本语法是ASCIIToRaster_conversion (in_ascii_file, out_raster, data_type)。我们可以在ArcGIS Desktop帮助文档里找到。

此外，gp还有一些专用的script命令，比如ListRasters、GetMessage等，大家可以在ArcGIS Desktop帮助文档查询Geoprocessor object。这些命令对于进行批处理迭代时十分有用。

在以上代码中，与geoprocessor有关的代码有，

gp=arcgisscripting.create()
gp.workspace=datadir
gp.ASCIIToRaster(txtfiles[i], grdnames[i], "FLOAT")

workspace设置了当前的工作目录，相当于在command window里执行 “workspace 目录”。

ASCIIToRaster执行实质的转换工作。如果此处我们需要执行位于ArcGIS扩展模块的命令，则需要在调用命令前，先对扩展模块（如spatial analysis）进行checkout。使用以下命令：

gp.CheckOutExtension("Spatial")
gp.CellStatistics_sa(...)

否则CellStatistic_sa会抱怨说没有License可用。

代码其余部分是常规Python代码。以下取得指定目录下的txt扩展名的全部文件名，并生成对应的grid文件名。

children=os.listdir(datadir)

#only ".txt" files
txtfiles=[]
grdnames=[]
for child in children:
    root,ext=os.path.splitext(child)
    if (ext=='.txt'):
        txtfiles.append(child);
        grdname='n'+root[8:19]
        grdname=grdname.replace('.','')
        grdnames.append(grdname)

Python十分有趣，如果有其它语言基础，学习起来并不难。对我们写ArcGIS脚本，有一般的Python知识就足够了。尽管效率不如传统的arcinfo AML高，但在功能上则不逊色。如果需要实现十分复杂的GIS功能，则建议使用c++, vb或.net语言操作ArcObjects。

When I install ArcGIS 9.2 developer SDK for .net framework 2.0 SP4, it will be stuck with time remaining 0 sec. So strange. It turns out to be caused by VC# express edition. A trick is to remove VCExpress process in Task Manager, thanks to contribution from ESRI forum posting. The installation will then continue.

Stupid ESRI.

I met a problem when I were running a routine in Command Line Window.
It closed the ArcMap. After then when I open ArcToolbox, everything gets empty.
See attached pic.

Many effort to figure it out by looking through Internet. Solution is rather simple.
Go to the register, navigate to ESRI fold under HKEY_CURRENT_USER/Software
then rename ESRI to ESRI_old. (Actually you can delete it safely.)

Now they get back.

Part III Advanced areal weighted interpolation

Zhuotong Nan (南卓铜) (zhn1@pitt.edu)

Motivation and objective

ArcView 9.2 comes with some basic interpolation approaches, for example, inverse distance weighted, Kringing, etc. In some cases, we want interpolate based on areal weighted. For example, Radar measurement can produce rainfall data over a study area. In order to use those Radar rainfall in hydrological model, e.g., VIC model, with a 1/8 degree spatial resolution, we need to do interpolation from a higher resolution, e.g., 4km by 4km. As shown in Figure 1, each grey cell has a measured value from Radar. We would like to generate an appropriate value for a VIC cell which might cover several Radar cells. In this part, we will illustrate an areal weighted approach.

Figure 1 Radar rainfall cells and 1/8 deg VIC cells

1. Create VIC grid polygons

Install Fishnet script with ArcView 9.2

ArcView 9.2 comes with a fishnet function (Data Management Tools → Feature Class → Create Fishnet) which is capable of generating a regular grid lines. However due to the limitation of ArcView license, it’s hard to implement conversion from polyline to polygon which is necessary to do areal weighted interpolation to be introduced in this section. Therefore, we will employ a free yet powerful script from ArcScripts website.

• Get this free script from http://arcscripts.esri.com/details.asp?dbid=12807. Unzip it to any location. We will see fishnet.dll in this location.
• Start ArcMap with a new document.
• Open Customize… dialog (Tools → Customize…), in Commands tab click Select from files… button to load fishnet.dll. You will see,

Figure 2 Customization dialog

Drag Create Fishnet to any place on toolbar. You will see on toolbar. Close Customize dialog.

Note, to uninstall the script, you can open Customize dialog again, and drag fishnet icon on toolbar to any place except toolbar.

Create 1/8 degree VIC grid polygons

VIC is a macroscale hydrologic model that solves full water and energy balances. It runs on basis of VIC grid. See http://ldas.gsfc.nasa.gov/LDAS8th/LDASspecs/LDASspecs.shtml for details of grid extent. We will generate VIC grid polygons using Fishnet utility download from ArcScripts in a scale of 1/8 degree.

• Open Create Fishnet utility ().
• Configure parameters as following figure.

Figure 3 Parameters for creating 1/8 degree VIC grid

• Wait its completion. This might take a couple of minutes. A polygon shapefile “vicgrd8.shp” in geographic coordinate system (i.e., longitude and latitude, specified by GCS_North_American_1927 parameter) will be generated upon its completion.

Reproject VIC grid polygons to Albers equal-area projection system

• Run Feature project utility under Data Management Tools → Projections and Transformations → Feature → Project in ArcToolbox. Configure parameters as following,

Figure 4 Reproject to Albers Equal Area projection system

The produced vicgrd8_alb.shp now is in Albers Equal Area projection system. Note here we specified a transformation method to convert NAD_1927 to NAD_1983. Remember vicgrd8.shp is in NAD_1927 system and now vicgrd8_alb.shp is in NAD_1983 system.

• Close ArcMap without saving the current document.

2. Generate polygons for XMRG radar data file

• Get example data from http://www.nws.noaa.gov/oh/hrl/dmip/2/data/0703200407z.asc.gz. Unzip it (using WinRar, etc) to our current workspace (c:\Workspace\Lab2).
• Convert the generated Grid in HRAP coordinate system to Stereogrpahic Polar coordinate system by using following equations. See this link http://www.nws.noaa.gov/oh/hrl/distmodel/hrap.htm for more details.

xster=hrapx*4762.5 - 401*4762.5
yster=hrapy*4762.5-1601*4762.5

Copy 0703200407z.asc to 0703200407z_ster.asc. Then open 0703200407z_ster.asc using any text editor, change its header to,

ncols 335
nrows 159
xllcorner -161925
yllcorner -6372225
cellsize 4762.5

Save 0703200407z_ster.asc.

• Convert 0703200407z_ster.asc to binary Grid file. In Command Line Window, run

ASCIIToRaster c:\workspace\lab2\0703200407z_ster.asc c:\workspace\lab2\g0703200407z FLOAT

• Define Polar Stereographic Coordinate System for g0703200407z grid file.

In ArcCatalog, right click g0703200407z to open its Properties window. Navigate to Spatial Reference section, click Edit… to define its spatial reference. Choose Define the coordinate system interactively. Select Stereographic (Polar view) and configure it as shown below.

Figure 5 Parameters for Polar Ster projection system

• Convert raster to polygon

The conversion from raster to polygon only support grid file in integer. To do this, run following commands in Command Line Window.

Times_sa c:\workspace\lab2\g0703200407z 100 c:\workspace\lab2\g07032004_100

Int_sa c:\workspace\lab2\g07032004_100 c:\workspace\lab2\g07032004int

At this point, the produced g07032004int is already an integer grid file.

RasterToPolygon c:\workspace\lab2\g07032004int c:\workspace\lab2\g0703200407_poly NO_SIMPLIFY VALUE

A field “GRIDCODE” in the generated polygon (g0703200407_poly) is corresponding to cell values in g07032004int, which is 100 times real rainfall value.

• Project XMRG polygons to Albers projection system

Using Data Management Tools → Projections and Transformations → Feature → Project in ArcToolbox to project g0703200407_poly in Polar Stereographic coordinate system to g0703200407_alb in Albers Equal Area Coordinate System. This target coordinate system can be imported from vicgrd8_alb.shp.

3. Calculate area contribution to each VIC cell

Calculate area for each VIC cell

Using CalculateAreas command to calculate area for each VIC cell in Command Line Window.

CalculateAreas c:\workspace\lab2\vicgrd8_alb.shp c:\workspace\lab2\vicgrd8_alb_area.shp

The field “F_AREA” in vicgrd8_alb_area.shp indicates area in square meter.

Overlay

Intersect utility is used to calculate areal contribution from XMRG polygons to each VIC cell.

• First we will clip the whole VIC grid shapefile with XMRG polygons extent. In Command Line Window,

Workspace c:\workspace\lab2
Clip_analysis vicgrd8_alb_area g0703200407_alb c:\workspace\lab2\vicgrd8_clp.shp

• In Arctoolbox, run Analysis Tools → Overlay → Intersect with input parameters specified as below.

Figure 6 Intersect input parameters

• Start ArcMap with a new document. Add intersect.shp. Save map document into c:\workspace\lab2\partIII.mxd.
• Open intersect.shp attribute table.

Figure 7 Attribute table of intersect.shp

Here FID_vicgrd and F_AREA come from VIC cells, while FID_g07032 and GRIDCODE come from XMRG polygon file. FID_vicgrd uniquely represents each VIC cell whose area is show in F_AREA field. GRIDCODE is actually equal to 100 * rainfall.

Querying attribute, for example FID_vicgrd = 1400, will extract all small polygons consisting a VIC cell.

Calculate areal weighted rainfall

• Add a new field VIC_AREA with values come from F_AREA.

Open intersect.shp’s attribute table, Click Options → Add Field…

Figure 8 Add a new field

Click the header of VIC_AREA column to select this field, right click it then select Field Calculator…, Yes upon warning message. Input [F_AREA] in the lower textbox then click OK to start calculation.

Figure 9 Field calculator

• Calculate area by running CalculateAreas in Command Line Window.

CalculateAreas intersect c:\workspace\lab2\inters_area

• Add a new field named “percentage” in double type to inters_area.shp.
• Calculate [percentage] = [F_AREA] * [GRIDCODE]/100 / [VIC_AREA].

Figure 10 Calculate Percentage field

• Export areal weighted rainfalls upon each VIC cell

Make inters_area layer visible, and make sure there is no selection in this layer. Otherwise use on the Tools toolbar to clear any existing selection. In ArcToolbox, run Analysis Tools → Statistics → Summary Statistics to calculate rainfalls cell by cell.

Figure 11 Summary statistics

Link interpolated rainfall data to VIC cells

• Make vicgrid8_clp.shp visible and be topmost. Link this layer’s attribute table to rainfalls.dbf by using context menu → Joins and Relates → Join. In Join Data window, enter inputs as shown below,

Figure 12 Join to rainfalls.dbf

• Export vicgrd8_clp to vicgrd8_rainfall.shp. The output shapefile will physically contain rainfall values (stored in field SUM_percen).
• Save partIII.mxd.
• Project vicgrd8_rainfall to rainfall_geo.shp.
• Add an alias “Rainfall” to field “SUM_Percen”. This makes the field name more straightforward.

Figure 13 Make field name self-explanatory

In order to use in combination with other raster, we can also convert this polygon shapefile to raster using FeatureToRaster command in Command Line Window.

4. Visualization

• Open partIII.mxd, switch all layers off except rainfall_geo. Set current data framework coordinate system to GCS_North_American_1983 using menu View→Data framework properties…→Coordinate System.
• Open rainfall_geo layer properties window, make symbology settings same as below. You can also choose a different color schema for this.

Figure 20 Symbology

• Export map for submission. See Figure 21 for reference.

Figure 21 Exported map

Conclusion

Although ArcView 9.2 has already included some basic interpolation approaches, it turns out in many cases hard to meet our specific requirements. This handout walks through how to interpolate Radar rainfall data to 1/8 degree VIC grid, illustrating necessary steps to implement areal weighted interpolation approach using ArcView. This handout also show us essential GIS operations such as project, clip, intersect, export, manipulate attribute table, link feature class to external table, etc, as well as making use of third party script.

(Http://nzt.spaces.live.com)

(转载请保持信息的完整，并加以适当引用)

ArcGIS command扩展：ArcGIS 9.2/C# 2的实现
Zhuotong Nan (南卓铜) (zhn1@pitt.edu)

(上节介绍了如何使用VS 2005 C# express来帮助我们生成ArcGIS Command扩展的代码框架。)

下面我们来介绍框架代码各部分的具体含义。

/// <summary>
/// Command that works in ArcMap/Map/PageLayout
/// </summary>
[Guid("eb9b3424-41d4-44e9-8ab0-c4740339db95")]
[ClassInterface(ClassInterfaceType.None)]
[ProgId("MyArcGISClassLibrary.BatchedChangeDatasource")]
public sealed class BatchedChangeDatasource : BaseCommand

ArcGIS是基于COM构建的，它只支持同样符合COM规范的动态链接库（Class library）。Guid必须不与任何一个COM组件相同。默认会生成一个新的全局标识符。我们也可以使用ArcGIS菜单里Developer Tools下面的Guid Generator来获取新的标识符（图10）。

图10，ArcGIS为我们提供了生成Guid的工具

ClassInterface和ProgId仍然是作为将.net类库导出成COM需要的一些属性。seal关键词指定了本命令不能再被继承。需要注意的是此类从BaseCommand上继承过来，BaseCommand实现了ICommand接口。9.2里BaseCommand放到ADF.BaseClasses命名空间下面，在之前版本是放在Utilities.BaseClasses下面。同时一些BaseCommand有关的辅助命令现在也放到了ADF.CATIDs空间。Utilities命名空间仍然存在，但当使用该空间下的BaseCommand在编译时会有deprecated的警告。

COM Registration Function(s)折叠起来的代码用于实现DLL到COM的自动注册和卸载过程。对于一般的应用，使用其生成代码就足够。

public BatchedChangeDatasource()构造函数里初始化了一些Command必要的信息，比如该Command所在的类别，这是以后在使用命令的时候存放的类别（如图11），同时也连接了一个命令表示的图标(m_bitmap)。这些变量是其基类BaseCommand里继承过来的。

string bitmapResourceName = GetType().Name + ".bmp";

如果不使用默认图标，则需要更改这里的名称。

Overriden Class Methods区域是需要开发者具体完成实现的地方。OnCreate在命令附加到ArcGIS时调用，一般这里传入的Object就是具体的Application，这里将Application对象存储到m_hookHelper.Hook属性里。通过Hook属性可以与ArcGIS Application进行必要的互访问。默认的OnCreate功能十分简单，只是判断了Application是否打开，Application是否有活动的View，如果有则使本命令可用。更多的判断则需要开发者来实现。

OnClick()是用户在点击Command后进行的操作。本案例中进行批量更改DataSource等功能便需要在这里实现。

此外，框架代码还将Project Build Properties设置成 Register for COM interop。对于手工生成的项目，必须记得设置此项。

图11，Project Build Properties设置

在具体添加代码前，我们先编译一下，看是否有问题。不幸的是默认的框架编译时会提示缺少引用。我们添加引用后（图13）编译通过。

图12，缺少引用错误

 
图13，增加ESRI.ArcGIS.Display引用

在类最前面添加私有变量声明。

public sealed class BatchedChangeDatasource : BaseCommand
{
    private ESRI.ArcGIS.Framework.IApplication _app;
    private string _layerNameToBeReplaced = "nwbi3";
   

_app接管Hook对象保管Application。_layerNameToBeReplaced指定要被替换Data Source的Layer名。同名的Layer会存在于不同的Data Frame里面。

因为我们需要实现的功能只对ArcMap有用，所以需要判断宿主应用程序是否MxApplication，如否，则禁用本命令。将OnCreate里的下面代码：

if (m_hookHelper == null)
        base.m_enabled = false;
else
        base.m_enabled = true;

更改成：

if (m_hookHelper == null)
    base.m_enabled = false;
else
{
    _app = hook as ESRI.ArcGIS.Framework.IApplication;
    if (_app is ESRI.ArcGIS.ArcMapUI.IMxApplication)
        base.m_enabled = true;
}

现在主要的功能是在OnClick里实现。要完成批量替代Data Source的功能，一般的思路是这样：在当前的Document（一个.Mxd文件即Document）里找到各个Data Frame（在ArcObjects，每个Data Frame就是一个Map），迭代各Map。在每个Map里迭代各Layer，判断此Layer是否具备指定的_layerNameToBeReplaced的名称，如是，将该Layer的Data Source代换成指定的。为了实现时间序列的代换，则还需要先形成目标Dataset Name的List，然后逐一代换。以下代码实现了将一个给定Layer Name的Feature层替换成固定的某个存储在File Geodatabase里的FeatureClass。targetGDBpath指定File Geodatabase的路径，具体的FeatureClass由 fcn.Name指定。如果此FeatureClass是存储在Geodatabase里的某个Feature Dataset里，可由fcn.FeatureDatasetName来指定。最终的Data Source的更换由IDataLayer.Connect来实现。

public override void OnClick()
{
    ESRI.ArcGIS.ArcMapUI.IMxDocument pDoc = (ESRI.ArcGIS.ArcMapUI.IMxDocument)_app.Document;
    if (pDoc == null) return;

    //_layerNameToBeReplaced ="nwbi3";
    string targetGDBpath = "c:\\temp\\fileGDB.gdb";

    ESRI.ArcGIS.Carto.IMaps maps = pDoc.Maps;
    ESRI.ArcGIS.Carto.IMap map;
    ESRI.ArcGIS.Carto.ILayer layer;
    ESRI.ArcGIS.Carto.IDataLayer dLayer;
    ESRI.ArcGIS.Geodatabase.FeatureClassNameClass fcn = new ESRI.ArcGIS.Geodatabase.FeatureClassNameClass();

    string conTofGDB = @"DATABASE= " + targetGDBpath;

    ESRI.ArcGIS.Geodatabase.IDataset ds = fGDBWorkspaceFromString(conTofGDB) as ESRI.ArcGIS.Geodatabase.IDataset;

    ESRI.ArcGIS.Geodatabase.IWorkspaceName wsName = ds.FullName as ESRI.ArcGIS.Geodatabase.IWorkspaceName;
    fcn.WorkspaceName = wsName;
    //IDatasetName dsName = ds.Subsets.Next().FullName as IDatasetName;
    //fcn.FeatureDatasetName = dsName;
    fcn.Name = "NWBI3";

    for (int i = 0; i < maps.Count; i++)
    {
        //each data frame
        map = maps.get_Item(i);
        for (int j = 0; j < map.LayerCount; j++)
        {
            layer = map.get_Layer(j);
            if (layer.Name == _layerNameToBeReplaced)
            {
                dLayer = layer as ESRI.ArcGIS.Carto.IDataLayer;
                if (dLayer != null)
                {

                    bool rev = dLayer.Connect(fcn);
                }
            }

        }
    }

    pDoc.ActiveView.Refresh();
}

COM风格的编程经常会看到很多的接口，不同的接口间不停的转换。其实接口一般是按功能进行分类，一个类通常实现了多个接口。比如FeatureLayer会实现不同的例如IDataLayer, IFeatureLayer, IDataset等，IDataLayer负责与数据源有关的操作，IFeatureLayer负责跟具备FeatureClass有关的操作，而IDataset实现了FeatureLayer作为Dataset有关的一些操作。

对Interface的熟悉应用建立在对相关类熟知的基础上。.net help十分全面，查找某接口，可以回溯到应用该接口的全面Class，而每个Class也会给出其实现的各个接口。通过这种方法，我们可以在一定程度上知道哪些接口是可以相互转换的。比如，IMap.get_Layer()查找帮助我们知道是返回ILayer，再结合我们期望是找到一个FeatureLayer，所以寻找FeatureLayerClass的实现的接口，发现其与数据源有关的操作封装在IDataLayer和IDataLayer2里。于是我们将layer cast成IDataLayer，再进行Connect操作。

以上代码中，fGDBWorkspaceFromString方法实现如下，作用是通过构建连接字符串来生成IWorkspace。

public ESRI.ArcGIS.Geodatabase.IWorkspace fGDBWorkspaceFromString(String connectionString)
       {
           ESRI.ArcGIS.Geodatabase.IWorkspaceFactory2 workspaceFactory;
           workspaceFactory = (ESRI.ArcGIS.Geodatabase.IWorkspaceFactory2)new ESRI.ArcGIS.DataSourcesGDB.FileGDBWorkspaceFactoryClass();
           return workspaceFactory.OpenFromString(connectionString, 0);
       }

至此，我们完成了一个很简单的可运行的Command扩展。回到构造函数，将m_category赋为“tong's Extensions”。编译后，打开ArcMap，打开tools菜单下的Customize...，弹出窗口如图14。左边的Categories可以找到tong's Extensions，右边的Commands可以看到我们的BatchedChangeDatasource。选定该命令，可以拖拉到任何工具栏。Close此Customize窗口后，命令就可以如内置一样运行。

 
图14，ArcMap的定制窗口

剩余的问题

对前面设定的期望目标，我们可以通过以下方法形成一个时间序列Dataset的名称。

private List<string> TargetRasterNames()
{
    //DateTime start = new DateTime(1997, 1, 27, 0, 0, 0);
    DateTime start = new DateTime(2002, 8, 19, 0, 0, 0);
    //DateTime end = new DateTime(1997, 1, 29, 0, 0, 0);
    DateTime end = new DateTime(2002, 8, 20, 0, 0, 0);
    List<string> rev = new List<string>();
    for (DateTime dt = start; dt < end; dt = dt.AddHours(1) )
    {
        //string dtStr = string.Format("c{0:MMddHH}c_3hm", dt);
        string dtStr = string.Format("{0:yyyyMMddHH}.asc", dt);
        rev.Add(dtStr);
    }
    return rev;

}

在Connect的时候，根据次序可以赋给不同的RasterName。在更换数据源时，如果确保都在同一Workspace下，则可以通过下面代码简单代换。

if (layer.Name==_layerNameToBeReplaced)
{
          dLayer= layer as IDataLayer;
          if (dLayer!=null )
          {
                 IRasterDatasetName pName = dLayer.DataSourceName as IRasterDatasetName;
                 IDatasetName dsName = pName as IDatasetName;
                 ii++;
                 dsName.Name = targets[ii];
            }
}

如果不在同一Workspace，IName还需要指定IWorkspaceName。

结论

本文结合实际工作的一个案例，演示了如何应用Visual Studio 2005 C# Express来开发扩展ArcGIS Command。以上演示还有很多需要改进的地方，比如要代换的Layer名，不能内置在代码里，我们可以通过生成一个Dialog让使用者来决定替换哪一层，再比如，我们可以设置更漂亮的图标以区别于已有的图标。但这些都已经不是本tutorial的内容。本tutorial的目的是希望大家通过阅读本文可以初步掌握如何使用C#进行ArcGIS扩展，大家更多需要关注的是OnClick里的Business logic。

有什么问题有与我联系或留言。转载请保留帖子全部信息。

ArcGIS command扩展：ArcGIS 9.2/C# 2的实现
Zhuotong Nan (南卓铜) (zhn1@pitt.edu)

为什么要做扩展?

很简单，ArcGIS提供了通用的功能，但不可能提供每个大家需要的功能。而且，对一些专业（如水文）也许最实用的功能都没有。这时候便需要对现有的ArcGIS进行扩展。

如何扩展?

有很多方法了，比如AML来扩展ArcInfo Workstation——我以前写过一个AML的简要指南，Google应该能找到。在Desktop下，可以用VBA扩展，或者用任何支持COM的现代语言来扩展。这里用Visual Studio C# 2005 Express来开发。对于独立应用程序，应当购买ArcEngine。ArcEngine可以认为是ArcObjects的一个子集。

准备工作

需要ArcGIS desktop .net developer kit。如果最初安装ArcGIS 9.2的时候没有安装上，那么到安装/删除程序里去重新修改安装。注意事先机器上得已经存在.net framework SDK，如果没有，安装ArcGIS看不到.net developer kit选项。.net SDK不同于平时大家安装的redistribution .net framework runtime。请到微软官方网站下载。

安装Visual Studio C# 2005 Express，建议在安装ArcGIS desktop .net developer kit之前进行，ArcGIS desktop .net developer kit将在Visual Studio C# 2005 Express上生成一些项目模板（如ArcMap class library），应用这些模板进行开发会节省很多时间。

我们的工作

如图1所示，要在一个Document里生成多个Map（多个Data Framework），Map里的内容是时间序列的小时降水量（不同的Raster Dataset）。右图仅演示了1天24小时的图。我们的挑战是要做10天，而且必须各个图采用同样的图例。

 
图1， 小时降水量图，我们希望实现的目标

解决方法

先设计好一个Data Frame，设置好Legend（不能是sketch，否则系列图就没有意义）。然后复制24份Data Frame，应用ArcMap的页面排版功能进行排版。然后将此文档复制10份（10天）。现在的问题是如何去更改每个Data Frame里的Data Source。我们需要的是实现一个批量处理的方法。采用ArcObjects进行二次开发，生成一个Command来完成这个事情，是个不错的方案。

建立一个新项目，点击File>New Project...（图2）。因为我们只想进行基于ArcMap的扩展，这里选择ArcMap Class Library（图3），敲入项目名称“MyArcGISClassLibrary”（不包括引号）。一般而言，项目按功能进行归类，一个可以实现很多的Command或者Tools扩展。所以没有必要对每个Command或者Tool都建项目。

图2，新建项目

 
图3，新建ArcMap Class Library

在随后出现的ArcGIS Project Wizard（图4），添加必要的引用。这里选择ESRI.ArcGIS.ArcMap，ESRI.ArcGIS.ArcMapUI，ESRI.ArcGIS.Carto，ESRI.ArcGIS.DataSourceGDB，ESRI.ArcGIS.DataSourceRaster，ESRI.ArcGIS.Framework，ESRI.ArcGIS.Geodatabase，ESRI.ArcGIS.System。这些都是常用的库，添加的这些库，会出现在VS项目的Reference里。如果这里不添加，也可以以后根据需要再行添加（如图4.1）。需要注意的是，这里用的ESRI.ArcGIS.DataSourceGDB是9.2新出来的File Geodatabase，据称比Personal Geodatabase性能要好，而且有更大的存储许可。点击Finish结束。

 
图4，ArcGIS Project Wizard

 
图4.1 可以从Solution的这个位置添加引用

该向导帮助我们添加了必要的Reference，并新建了Class1.cs。但我们不需要它。如图5所示，删除。

 
图5，删除Class1.cs

移动Solution Explorer窗口的项目上，如图6，添加一个新项，点击New Item...。在新弹出的窗口中选择Base Command（图7）。并设置名称为"BatchedChangeDatasource.cs"。在New Item Wizard Options里选择ArcMap Command（如图8）。

 
图6，新建Command类

 
图7，选择Base Command

 
图8，选择ArcMap Command

点击OK，Wizard生成足够多的东西（如图9），我们在此基础上再写一些代码以完成我们的预期功能。

 
图9，Wizard生成的代码框架按Ctrl+S保存此Solution。Solution名称设成MyArcGISClassLibrary.sln。

(下一节我们将介绍代码的作用，和如何添加自己的功能以实现既定目标。)

1. 关闭fgdc和geo-network同步器。在metata菜单里，使用set working synchronizer命令关闭两者。
2. 建立一个vba宏，删除xml里的fgdc标签。edn有现成可用的宏，在esri的edn开发网站上搜索 delete FGDC metadata，估计就可以找到。（或向我索要）建立好了，在ac里选择上需要处理的文档，支持宏即好。其实分析一下，就只有几个标签，手工去删除也不难。而且最好能有个自己的工具去批量删除。
3. 关闭 arcims server上的基于内容的 validation功能，用处不大，那是为metadata explorer设计的，如果不用explorer，这个功能就没多大用处了，而且此功能实质是检查 fgdc标签的内容（或者说iso的内容处在低于fgdc内容的优先级别上）。