Category Archives: GIS

MODIS过境时间及相关

Zhuotong Nan(南卓铜, [email protected])

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

na2002_08_02_214 na2005_12_08_342_aqua
Aug 2, 2002(左),Dec 8, 2005(右)

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

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

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

Kappa tool

Zhuotong Nan (南卓铜,[email protected])

网上可以找到的计算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 data1rldas%%i data2rldas%%i  -k:.kappa.txt
)
echo Done.

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

[]表示必须,{}表示可选。对ESRI grid,图像指定到目录名;对个人geodatabase(通过Access数据库保存)的Raster dataset,以geodb.mdbraster指定;文件geodatabase的Raster dataset,以geodb.gdbraster指定。如果指定日志文件,在屏幕上输出的内容将同时保存到日志文件中;如果指定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联系我。

在ArcGIS中应用脚本实现批处理

Zhuotong Nan(南卓铜) [email protected]

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 Filesarcgisbin目录下的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。

Stuck when updating ArcGIS 9.2 sdk for .net sp4

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.

Every tool is empty in ArcToolbox 9.2

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.

arctoolbox_problem

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.

太糟糕了

ArcGIS 9.0与IE 7冲突,表现在运行ArcToolbox的任何命令,都会关闭整个app。
以为是我机器的问题,figure it out花了好几个小时
后来发现是本身的问题,ArcGIS 9.1提供了Patch,ArcGIS 9.2已经兼容
但ArcGIS 9.0是manual support,建议将IE 7卸掉。
ESRI不是一般的黑。变着法子让大家花钱去升级。

See here for official statement.

Advanced Areal Weighted Interpolation by ArcView 9.2 – draft version

Part III Advanced areal weighted interpolation

Zhuotong Nan (南卓铜) ([email protected])

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.

clip_image002[8] clip_image004[8]

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,

clip_image006[6]

Figure 2 Customization dialog

Drag Create Fishnet to any place on toolbar. You will see clip_image008[14] on toolbar. Close Customize dialog.

Note, to uninstall the script, you can open Customize dialog again, and drag fishnet icon clip_image008[15] 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 (clip_image008[16]).
  • Configure parameters as following figure.

clip_image010[6]

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,

clip_image012[6]

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

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:workspacelab2703200407z_ster.asc c:workspacelab2g0703200407z 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.

clip_image014[4] clip_image016[4]

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:workspacelab2g0703200407z 100 c:workspacelab2g07032004_100

Int_sa c:workspacelab2g07032004_100 c:workspacelab2g07032004int

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

RasterToPolygon c:workspacelab2g07032004int c:workspacelab2g0703200407_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:workspacelab2vicgrd8_alb.shp c:workspacelab2vicgrd8_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:workspacelab2
Clip_analysis vicgrd8_alb_area g0703200407_alb c:workspacelab2vicgrd8_clp.shp

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

clip_image018[4]

Figure 6 Intersect input parameters

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

clip_image020[4]

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…

clip_image022[4]

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.

clip_image024[4]

Figure 9 Field calculator

  • Calculate area by running CalculateAreas in Command Line Window.

CalculateAreas intersect c:workspacelab2inters_area

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

clip_image026[4]

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 clip_image028[4] on the Tools toolbar to clear any existing selection. In ArcToolbox, run Analysis Tools → Statistics → Summary Statistics to calculate rainfalls cell by cell.

clip_image030[4]

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,

clip_image032[4]

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.

clip_image034[4]

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.

clip_image036[4]

Figure 20 Symbology

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

clip_image038[4]

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)

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

SimCity开源了

On January 10 2008 the SimCity source code was released under the free
software GPL 3 license under the name Micropolis.

http://en.wikipedia.org/wiki/SimCity (国内可能要加代理才能访问)
http://simcitysocieties.ea.com/about.php?languageCode=1

Extending ArcGIS Command: A case tutorial (II)

ArcGIS command扩展:ArcGIS 9.2/C# 2的实现
Zhuotong Nan (南卓铜) ([email protected])

(上节介绍了如何使用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)。

image
图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。对于手工生成的项目,必须记得设置此项。

image
图11,Project Build Properties设置

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

image
图12,缺少引用错误

image 
图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:tempfileGDB.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窗口后,命令就可以如内置一样运行。

 image
图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。

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

Extending ArcGIS command: a case tutorial

ArcGIS command扩展:ArcGIS 9.2/C# 2的实现
Zhuotong Nan (南卓铜) ([email protected])

为什么要做扩展?

很简单,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都建项目。

image
图2,新建项目

 image
图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结束。

image 
图4,ArcGIS Project Wizard

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

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

 image
图5,删除Class1.cs

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

image 
图6,新建Command类

image 
图7,选择Base Command

image 
图8,选择ArcMap Command

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

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

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