Tag Archives: tutorial

Matlab科学数据处理:如何补上缺失的时间序列

Matlab科学数据处理:如何补上缺失的时间序列

Zhuotong Nan (南卓铜, giscn@msn.com)

问题提出

假设我们有一个时间序列文件,内容如下

286880 19480223 -5.7

286880 19480227 10.8

286880 19480229 7.7

286880 19480301 -15.7

286880 19480302 3.2

286880 19480303 7.7

286880 19480305 21.2

第一列是站点名,第二列是时间信息,第三列是温度值。时间并不连接,我们经常需要把时间补全,这些时间上的观测以缺失值(比如999.9)替代。即,我们希望是,

286880 19480223 -5.7

286880 19480224 999.9

286880 19480225 999.9

286880 19480226 999.9

286880 19480227 10.8

286880 19480228 999.9

286880 19480229 7.7

286880 19480301 -15.7

286880 19480302 3.2

286880 19480303 7.7

286880 19480304 999.9

286880 19480305 21.2

那么在matlab里要如何实现呢。

基本思路 Continue reading

多元回归 / 多元拟合matlab实现

MatLab内置了很多拟合的功能,比如分布、指数、多项式等拟合。对应的函数是fit。基本用法是:

cfun = fit(xdata, ydata, libname)

其中xdata是自变量x序列,ydata是因变量,libname 是模型名,比如’power1’是 y= a* x ^b 形式,Matlab内置了很多模型,可以从帮助文档里看到。对于这些拟合,MatLab有个图形工作叫 cftool,在命令行里敲入可以打开。

然而我们有时候需要做二元、多元拟合。

二元拟合,即使用模型 z = f(x, y),在MatLab也叫面拟合 (Surface fitting),有专门的图形工具 sftool,命令敲入sftool可以打开。

但如果是多元拟合,这时我们需要懂一点代码的工作。

其实不管一元、二元,代码级别是一致的。使用的都是fit命令。

对于使用 y = a* x^b这样内置的指数一元拟合,我们可以看到其libname 叫’power1’。所以采用以下命令:

cfun = fit(x, y, ‘power1’)

即从已知x、序列是得到拟合函数。coeffvalues(cfun)可以得到拟合参数,即这个例子里的a、b。如果还需要得到拟合指标,使用[cfun gof] = fit(x,y, ‘power1’),返回的gof包括了R squre、RMSE等多种指标。

对于二元拟合,比如 使用 z= a + b*sin(m*pi*x*y) + c*exp(-(w*y)^2模型进行拟合,采用类似的代码:

ft = fittype( ‘a + b*sin(m*pi*x*y) + c*exp(-(w*y)^2)’, ‘indep’, {‘x’, ‘y’}, ‘depend’, ‘z’ );

opts = fitoptions( ft );

opts.Display = ‘Off’;

opts.Lower = [-Inf -Inf -Inf -Inf -Inf];

opts.StartPoint = [0.741864193163802 0.681571628918889 0.811658042300087 0.275716137656952 0.738979809065491];

opts.Upper = [Inf Inf Inf Inf Inf];

opts.Weights = zeros(1,0);

[fitresult, gof] = fit( [x, y], z, ft, opts );

其中fittype自定义了需要的二元函数。其中’indep’指明自变量、’depend’指明因变量。即通过fittype构造一个二元拟合模型,说明其中的自变量和因变量,定义拟合的选项,这些选项经常可以不要专门指定,采用默认的即可。然后使用 fit的扩展形式得到拟合函数。

对于多元拟合,一个函数是 regress多元线性回归,另一个是nlinfit 非线性回归。两者都在statistic toolbox里。regress的基本形式是

b = regress(y, X)

其中X是多个自变量,y是因变量。

对于多元非线性拟合,可以使用nlinfit,基本形式是:

beta = nlinfit(X,Y,modelfun,beta0)

其中modelfun构造非线性多元模型,比如以下例子:

beta = nlinfit(X,y,@hougen,beta0)

使用一个hougen函数作为拟合模型,这个函数可以自行构建。

Matlab也提供了非线性多元拟合的图形化工具nlintool。

此外类似的还有 fitnlm函数提供了类似nlinfit的非线性回归,用法与nlinfit类似,但对返回结果进行包装,可以更简便的使用诸如predict等函数对非线性模型进行预报计算。

Image(1)

制作符合学生信息表要求的数码照片

小学生最近填学生信息表,要求提交一张数码照片。朵朵小朋友的事情就落到我头上了。照片的要求还多多。

比如,1. 淡蓝色背景。2. 照片规格:358像素(宽)x 441像素(高),32mm x 26mm, 分辨率300dip以上。(原文如此)。3. jpg格式。4. 文件大小小于60K。

大家注意的32mm x 26mm没有,这里32mm和26分别是什么,都没有讲,与像素高宽对比,才能知道,32mm只能是高,26mm是宽。然后细心的可能会注意到另一个很脑残的事情,按358像素宽和26mm宽的比率,这个分辨率并非300dpi或者是一个整数,而是349.758xxxx。如果在后期处理的时候,不注意到这个分辨率的调整,是怎么也不可能同时满足像素和实际尺寸的要求的。脑残的是也不知道定一个整数的dpi。

下面讲如何在photoshop里进行简单的处理以满足这些要求。

1. 我的原文件是朵朵的护照照片,原大小是 602像素x602像素,300dpi。白底。

image

2. 打开原照片。在Image > Image size里进行调整。

确认constrain proportions打勾,在resolution里填那个脑残的比例,即349.758 pixel/inch(即dpi)。在document size的height里填 3.2 单位是cm。这时由于constrain proportions打勾,会按比例缩放。点击确定。

image

3. 在view > rulers里把rulers打勾,即显示标尺。选择剪裁工具,从最左边开始,拉选 358个宽的像素。如图。

image

不是严格的358没有关系,以后可以细调。按向右的方向键把框框移到合适的中间位置。如图。

image3

回车键确认。正式进行剪裁。

4. 重新打开Image > Image size,这时如图

image

注意到Width并非 358。把Constrain Proportions的勾去掉,把Width 改为 358。这时注意到 Document size里的Width 也相应变为 2.6 cm。确认即可。现在已经在正确的尺寸了。

5. 从网上找一个淡蓝的证件照。比如这个http://img.club.pchome.net/upload/club/2004/7/15/chenqingfa_1089898576.jpg。下载保存。在photoshop里打开这个证件照。用吸管取其背景淡蓝色。或者直接用 # 66ccff这个颜色值设置前景色。

6. 在朵朵小朋友的照片里,用魔术棒点了一下白色的区域,这时整个白色区域就全部被选上,如图。

image5

7. Edit > Fill 里,使用前景色填充,确认后如图。

image

8. 文件 > 另存为,取合适的文件名,在jpg 选项里,注意文件大小,拖拉 质量,使得右边上的文件大小小于60KB(可以适当大于60KB,保存后的实际大小会小于右边这个数字)。然后保存。

image

9. 恭喜,一个满足全部要求的照片完成了。

10. 不要忘记,学校是要求用身份证号来命名的。

ThinkPad x220 更换硬盘及数据迁移

Zhuotong Nan (giscn@msn.com)

原x220 4286-CTO 是320GB 7mm硬盘。已经没有空间了。买了一块500GB的日立 Hitachi Z7K500 (7200转)7mm硬盘。我此前写过一个x200更换硬盘并迁移数据的帖子。那里面描述的方法肯定是可以用的。但这次尝试另一种更为简便的方法。

1. 笔记本里原硬盘的数据必要备份。以免万一。

2. 拆硬盘的视频可以看这里,从2分05秒开始。准备后十字螺丝刀。

IMG-20120815-00009

3. 用Orico 6628 series 的tool free dual bay hard drive dock 进行盘对盘clone。注意,源盘和目标盘一定不要弄错了,所以在切换为clone模式,按下start前务请确认没有错。

IMG-20120815-00007IMG-20120815-00008
*一定注意源盘(source)和目标盘(target),不要放错。源盘要放原笔记本里的盘,目标盘是新的500G盘。

4. 开始clone,指示灯为闪红,大约需要30分钟。结束的时候有三声鸣音,指示灯变蓝后为合适了。

5. 把新盘装进x220。重启。

6. 进disk management,把未分配的空间并入原盘(my option),或者建新盘。如果built-in功能不好用,建议用acronis的相关工具。

7. 扩展后的结果如下图。

Image(4)

一切OK。

AHK实现批量网页数据提取

Author: Zhuotong Nan (giscn@msn.com)

AHK是很好用的在windows下实现自动化的脚本语言,对于有一定编程基础的童鞋,很快就可以上手写代码。ahk的官方网址是http://www.autohotkey.com/,其forum里有大量的他人写好的脚本。

本文实现一个从网页里自动提取数据的一个例子。

目的

从40407.com网站里自动批量得到《武林3》70区豪侠礼包代码。原始网站界面如下,

image

1. 进入http://fahao.40407.com/11417.html,点击“淘号”按钮,如果未登陆,则提示登录。登录后弹出淘号界面如下图,

image

从界面中获取以XX打头的15位代码。

2. 关闭“查看卡号”,然后重复步骤1,可以得到一系列的代码。

现在的目的是实现其自动化,并把得到的代码保存到一个文本文件(gifts.txt)里。

实现

1. 点击“我要淘号”时,本地与远程有一个通讯。利用chrome自带的inspect element的network功能(也可以用任何网络分析工具),在“我要淘号”上右键选 inspect element, 到Network tab上,这时点一下“我要淘号”执行一次通讯,Network将捕捉到通讯信息。

image

我们可以看到通讯的header。我们需要做的是模拟chrome的通讯用ahk来实现。

2. ahk有几个有用的他人写的http通讯脚本,一个叫httpquery,一个叫httprequest,都很优秀,但前者2010年后已经不再继续支持,并且不支持cookies文件头。而cookie信息是我们这个任务里的关键,否则“我要淘号”返回是要求登陆的信息。httprequest可以从这个网址下载:http://www.autohotkey.com/community/viewtopic.php?t=73040

3. 代码实现。本代码实现了自动获取500个号码。

; Author: Zhuotong Nan (giscn@msn.com)
; Date: 8/8/2012

loop, 500
{
    url := "http://fahao.40407.com/app/common.php"
    data := "a=taohao&id=11417"

    header :="Cookie: PHPSESSID=da317c26545caafa204520e4c89276c9; DedeUserID=427080; DedeUserID__ckMd5=d770e5fdd4d5a237; DedeLoginTime=1344523159; DedeLoginTime__ckMd5=2233e446322305d5; Hm_lvt_75e3d45d3ad80ecb1af209e743d1b514=1344521115290,1344521187845,1344523168489; Hm_lpvt_75e3d45d3ad80ecb1af209e743d1b514=1344523168489`nReferer: http://fahao.40407.com/11417.html`nConnection: keep-alive"
    options :=""
   
    bytes := httprequest(url, data, header, options)

    sleep 500

    ;parse gift code
    pos:=Regexmatch(data,"XXw{13}")

    if pos!=0
    {
        code := substr(data,pos,15)
        fileappend,%code%`n,gifts.txt
    }
   
    sleep 1000
   
    fileappend,%data%`n,test.html
}

#include httprequest.ahk

ESC::
exitapp

 

4. 注意httprequest有4个参数,参数1是访问的地址,参数2是上传的数据(post data)和下载过来的内容(放同一个参数里),参数3是上传的http头,响应的http头也在参数3里,参数4是选项。我们需要使用的是前3个参数。注意参数3 header的构造,不同的key-value对使用分行符(`n)来分隔。而且这里必须要增加cookie信息,这个信息可以从chrome里分析得到。传回来的数据存放在参数2,data里。

5. 使用了正则表达式,从下载过来的data 里进行pattern分析。礼包代码是以XX开头的15长度的大小写字母和数字组成。表达式为XXw{13}。

6. #include httprequest.ahk 把httprequest脚本包括进来,所以前文可以调用httprequest函数。

7. 获取的代码存放在 gifts.txt里。

image

8. 使用本代码前,必须事先用chrome (或者IE等别的浏览器,需要对应获取其cookie)登录成功。按esc键随时退出。

结论

本文演示了使用AHK从网页自动批量获取指定数据。使用了httprequest方法,其中关键是必须维持cookie信息。

A Lab for using the PIHMgis and PIHM model (II)

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.

clip_image002

Figure 3 QGIS Plugin Manager

Now the PIHMgis toolbar appears as shown in Figure 4.

clip_image004

Figure 4 the PIHMgis v2.0.0 toolbar

3. Load DEM to PIHMgis for displaying.

Layer > Add a Raster Layer…, select lab3.shale.hillsDEMsh_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),

clip_image006

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).

clip_image008

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.

clip_image010

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.

clip_image012

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 clip_image013 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.

clip_image015

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).

clip_image017

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.

clip_image019clip_image021clip_image023

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.

clip_image025

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.

clip_image027

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.

clip_image029

Figure 14 Catchment grid generation interface

clip_image031

Figure 15 Convert catchment grid to vector

clip_image033

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.

clip_image035

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.

clip_image037

Figure 18 The Attributes toolbar

clip_image038

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.

clip_image040

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.

clip_image042

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.

clip_image044

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.

clip_image046

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.

clip_image048

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)

clip_image050

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.

clip_image052

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

clip_image054 clip_image056

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

Close QGIS 1.0.

A Lab for using the PIHMgis and PIHM model (I)

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.

clip_image002

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.0bin”, 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.

clip_image004

Figure 2 the QGIS 1.0 interface

使用matlab拟合Gamma分布

南卓铜 (zhn1@pitt.edu)

任务描述: 有一批数据,需要将其histogram拟合成Gamma分布,histogram分20个bin,从0到100。结果如下图所示:

hist_dist_final

实现步骤:

1. 装载数据进来,数据文件是data.for.hist.txt,文件结构是一列无header的数据,共43个数据。

>> data=load(‘data.for.hist.txt’)

data是43×1的矢量。

2. 打开distribution fitting tool

>> dfittool

选择Display type为Density(PDF)。点击Data…,弹出Data对话框。

在此对话框内,设置Data为 ‘data’矢量,Censoring和Frequency为 none。设置data set name为dataset,然后点击 Create Data Set。

点击Set Bin Rules,设置Bin width为5。其余不变,OK确认关闭Set Bin Rules。这时,我们将data的histogram设置为每个bin长度为5。在不同的应用里,bin根据具体情况调整。关闭Data对话框。

3. Gamma拟合

点击New Fit…,在New Fit对话框里选择Data为dataset,Distribution为Gamma,点击Apply进行拟合。点击Close关闭New Fit对话框。

此时,看到下图的效果:

image

我们需要修改它的X域,使之匹配我们的histogram的X最大值100。点击Tools菜单下的Axes Limit Control,将X Lower Limit和X Upper Limit更改为 0和100。

image

4. 导出成figure

使用File菜单下的Print to Figure导出成figure。对figure进行必要修改。

>>figure(1)  (注,如果已经打开多个figure窗口,导出的figure序号可能不是1而是其它数字,注意figure窗口标题上的Figure n里的n)
>>legend hide
>>xlabel ‘1/32 degree, mm/day’
>>ylabel ‘Frequency’

此时效果如下,

image

5. 修改figure格式

设置figure为Tools-> edit plot状态,打开View->Property Editor。

将histogram的plot type从line改成Area,配置合适的face color和edge color。

将x轴的范围X Limits设置为0 to 100。

image

这时,histogram压着拟合曲线,选择拟合曲线,右键Cut后,再Paste进来,可以将拟合曲线带到最上面。

6. 修改数据

注意到目前为止,Y轴仍是Density。我们需要将Y坐标修改为Frequency(或者count)。Count=Density x bin width x data count,这里bin width是5,data count是43。因此,我们需要对Ydata乘以5×43=215。

选择拟合曲线,

>>ydata=get(gco,’YData’)
>>ydata_215 = ydata*215
>>set(gco,’YDATA’,ydata_215)

选择histogram,

>>set(gco,’YDATA’,get(gco,’YDATA’)*215)

修改xlim和ylim

>>ylim([0 10])
>>xlim([0 100])

将figure还原为非编辑状态。Tools->Edit Plot,并关闭Property editor。

image

以上步骤中数次修改xlim和ylim,并不是每次都需要,这里仅是为了演示需要。

7. 保存图形

>>print -dpng ‘histogram_dist.png’

在figure对话框的File-> Save As…保存为MATLAB Figure(.fig),以后可以通过File->Open…打开。

结论

本文演示了如何使用Matlab对数据的histogram使用Gamma分布进行拟合。介绍了matlab的常用命令,Distribution Fitting Tool,figure的格式化编辑功能,对figure对像数据的修改等。也介绍了比如如何调整figure里对像的z-order等实用技巧。

http://cid-0ea641a5a7f665a1.skydrive.live.com/embedrow.aspx/Public/data.for.hist.txt

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

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

Extending ArcGIS Command: A case tutorial (II)

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)。

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。

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