GIS中公交数据的下载和处理

　　前两天讲座的Demo中，演示了一个可以在等公交时查看将要到站车辆信息的移动应用，其中用到了北京若干条公交线路的数据。不少朋友都想知道这些数据是如何获取的，我把方法共享出来，希望能为大家GIS应用的道路上扫清一些不必要的障碍。
　　为了获取公交数据，我对比了几个具有查询公交线路功能的在线地图网站，查看http的请求结果，发现百度地图的http响应最接近成品状态，它的结果是标准的json数据，含有查询公交线路上所有的站点信息和公交路线信息等。简单验证了一下，如果从百度地图网站上发起请求，可得到所有明文信息，如果从其他站点发起请求，所有涉及到地理坐标的值会被加密。
　　我想完成的事情是，获取关心的公交线路数据，然后把他们做成Geodatabase中的FeatureClass来使用。只要能得到有效的JSON数据，就可以利用ArcGIS中的Python脚本来完成我想要的结果。
　　第一步需要获取JSON数据。为了获取批量的结果，我们需要使用百度地图的API。我学习了酸奶小妹妹写的一篇现成的文章，知道了从Javascript中可获取我想要的JSON信息。下一步就是如何将js中的JSON结果保存成文本文件。js直接写txt不太好办，借助Silverlight与html内容的交互功能，就可以办到这点。我写了一个Silverlight应用来完成这个工作，首先填写想要查询的目标城市，比如“西安”，然后填写想要查询的公交线路，多条线路间需要用英文半角逗号隔开，点击“开始下载”就会从百度得到JSON结果，保存成txt文件以备后用。
　　第二步需要用Python处理JSON数据以生成FeatureClass。我写了一个python脚本来完成这个工作，并把它制作成了Script Tool，这样就可以通过调用GP工具的形式来生成我们的地理数据了。
　　第三步数据后处理。百度API中得到的地理数据是经纬度坐标，生成FeatureClass后可利用Project工具投影到我们想要的坐标系中，比如WGS 1984 Web Mercator。对于火星坐标，一个城市范围内的偏移可以认为是线性的，只需在Arcmap中参照底图整体平移即可。
　　比如，我在这个网站上参照页面源码可以得到西安市所有公交线路列表，简单处理后变成这样：

地铁1号线,地铁2号线,地铁3号线,地铁4号线,地铁5号线,地铁6号线,机场快轨,环山1号线,环山2号线,二环1号线,二环2号线,高新1号线,高新草堂专线,教育专线,五龙专线,泾渭环线,1,2,4,K5,6,7,K8,9,10,11,12,13,14/K14,15,16,17/K17,K18,19,20,20区间,21,22,23,24,25,26,27,28,29/K29,30/K30,31,32,33,34,K35,36,37,38,39,40,41,42,K43,44,45,46/K46,47,48,50,102,103,104,105,106,107,108,117,118,161,162,163,201,202/K202,K203,204,K205,206,207,207区间,K208,209,K210,211,212,213,214,215/K215,215区间,216,217,218,219,220,22,222,223,224,225,226,226区间,227,K228,229,230,230区间,231,232,233,234,K235,236,237,238,239,240,241,242,243,251,252,260,260区间,261,262,263,264,265,266,300/K300,301,302,303,307,308,309,311,312,313,315,316,318(车城),318(全程),319,321,322,323,324,规划325,规划326,328,330,331,336,351,规划356,K400,401,402,403,405,406,407,408,409,410,411,416,500,501,502,503,504,506,507,508,509,510,511,512,规划513,规划515,517,518,519,521,522,523,525,526,527,528,530,600,601,602,603,604,604区间,K605,606,606区间,607,608,K609,611,612,K616,K618,619,K630,630区间,K631,K700,701,702,703,704,705,706,707,708,709,710,711,712,713,714,715,716,717,718,719,720,721,722,723,728,729,800/K800,801,规划802,规划821,规划822,900,900区间,901,901区间,902,903,904,905,906,907,908,910,911,912,913,914,915,916,917,918,919,920(汤峪),920(焦岱),921,922,923,924,926,927,928,929,930,东南专线,宏华线,安运线,西高客运,规划西阎线,通宵1号线,通宵2号线,通宵3号线,通宵4号线,游1,游2,游3,游4,游5(306),游6,游7,游8(610),游9,游10,盛唐之旅,曲江新景线,曲江旅游观光轻轨,世园环线,世园1号线,世园2号线,世园3号线,世园4号线,世园5号线,世园6号线,世园7号线,世园8号线,世园9号线,世园10号线,世园12号线,世园13号线,3-18,4-04,4-05,4-06,4-07,4-08,4-09,4-10,4-11,4-12,4-13,4-14,4-15,4-16,4-17,4-18,4-19,4-20,4-21,4-22,4-23,4-25,临潼101,临潼102,临潼201,临潼301,临潼826,阎良601,阎良901,高陵101,高陵102,高陵103,高陵105,高陵106,高陵108,高陵109,高陵201,高陵801,高陵802,高陵803,户县801,户县802,户县803,户县804,户县805,户县806,户县808,户县809,户县810,户县811,户县812,户县813,户县819,户县820,咸阳13,咸阳15,咸阳21,咸阳22,咸阳28,咸阳29,咸阳59,咸阳副59,

　　利用这个Silverlight应用即可下载到JSON结果，再用一下GP工具，就可生成西安市公交线路的地理数据。


　　当然你也可以把它导出成shapefile，方便其他同事，同学或老板使用。
　　GP工具，Python脚本，示例FeatureClass可在这里下载：http://arcgis.com/home/item.html?id=e0f8316d91fb43d49a81a76946f9a03c

Python魅力之地理数据快速处理方法

        Python作为一种成熟的脚本语言，发展势头迅猛，长期徘徊在“编程语言谱”排行前5的位置，在脚本语言中仅次于php。通过Python，开发人员可以快速试验一个想法、研究一种对象属性或不同算法，而无需编译、执行和测试任何代码。正是因为跨平台、简洁、优美的特点，它也如同GIS应用渗透于各各个行业一样，渗透在科学计算的各个领域。在GIS领域，几乎可以说开源言必称Python；而ArcGIS则将Python称为the scientific programming language for GIS：ArcGIS Server的猛将Geoprocessing Service依赖于Python，从9.4将Python IDE集成到ArcMap中也可见一斑。
        前一阵Flyingis贴出了一篇类似的文章，不过还是通过动手完成一个实例，来看看它在ArcMap中的数据快速处理能力吧。描述：上传到Panoramio的照片都会有位置信息（经纬度坐标），有些朋友非常可敬，上传的照片数量可达上万张。分析一下某个用户拍照的活动范围，是件有趣和有意义的事情。思路：1、利用Panoramio的API，获取指定用户的所有照片信息；2、将关心的信息，比如作者、照片名称、照片文件链接提取出来，并将其存储到地理数据库中。
        第一步：获取照片信息。先看一下Panoramio的API，REST风格，返回JSON字符串。比如通过http://www.panoramio.com/map/get_panoramas.php?order=popularity&set=public&from=0&to=20&minx=-180&miny=-90&maxx=180&maxy=90&size=medium，即可按点击次数排序获得世界范围内所有受欢迎的照片。

{
"count": 773840,"photos": [
{
"photo_id": 532693,
"photo_title": "Wheatfield in afternoon light",
"photo_url": "http://www.panoramio.com/photo/532693",
"photo_file_url": "http://static2.bareka.com/photos/medium/532693.jpg",
"longitude": 11.280727,
"latitude": 59.643198,
"width": 500,
"height": 333,
"upload_date": "22 January 2007",
"owner_id": 39160,
"owner_name": "Snemann",
"owner_url": "http://www.panoramio.com/user/39160",
},
{
"photo_id": 505229,
"photo_title": "Etangs près de Dijon",
"photo_url": "http://www.panoramio.com/photo/505229",
"photo_file_url": "http://static2.bareka.com/photos/medium/505229.jpg",
"longitude": 5.168552,
"latitude": 47.312642,
"width": 350,
"height": 500,
"upload_date": "20 January 2007",
"owner_id": 78506,
"owner_name": "Philippe Stoop",
"owner_url": "http://www.panoramio.com/user/78506"
}, ...
]
}


        Panoramio文档说有773840张，但是JSON里的count只有338张，搜其论坛，发现这是一个bug。当然，为了安全起见，一次请求最多只能返回100张照片的信息，Python刚好来完成若干个“100”的拼装任务。代码如下：

import arcgisscripting,json,urllib,os,threading
lphotos=[]#list saves all photo json objects
lock=threading.Lock()

#function to construct the panoramio request url and retrieve the corresponding photos as json array
#for eg: http://www.panoramio.com/map/get_panoramas.php?order=upload_date&set=139922&from=0&to=20&minx=108.355671725&miny=33.807692976&maxx=109.475020647&maxy=34.732003828&size=medium
#set: public (popular photos); full (all photos); user ID number
#size: original; medium (default value) ;small ; thumbnail ; square ; mini_square
def retrievephotos(set,fromindex,toindex,size):
global lphotos,gp
#extent=r"&minx=108.355671725&miny=33.807692976&maxx=109.475020647&maxy=34.732003828"
extent=r"&minx=-180&miny=-90&maxx=180&maxy=90"
url=r"http://www.panoramio.com/map/get_panoramas.php?order=upload_date&set="+set+"&from="+str(fromindex)+"&to="+str(toindex)+extent+"&size="+size
try:
jsonresult=json.read(urllib.urlopen(url).read())
except:
raise Exception,"get panoramio result error"
photos=jsonresult["photos"]
if len(photos)<toindex-fromindex:
print "warning: can't get photos from server "+str(fromindex)+"--"+str(toindex)
#gp.addwarning("warning: can't get photos from server "+str(fromindex)+"--"+str(toindex))#typeerror:'no type' object is not callable
else:
for photo in photos:
lock.acquire()
lphotos.append(photo)
lock.release()
return jsonresult["count"]#the total of photos of the specified set

#use multithread to retrieve photos, the number of threads=(end-start)/100
#start: 0 in usual
#end: the total of photos
def getphotolist(set,start,end,size):
list=[start]
while list[-1]<end:
list.append(list[-1]+100)
threads=[]
gp.addmessage(str(len(list))+" threads has launched...")
for i in list:
threads.append(threading.Thread(target=retrievephotos,args=(set,i,i+100,size),name=str(i)))
for t in threads:
t.start()
for t in threads:
t.join()#suspend the parent(main) thread until this thread is over

        在运行过程中，发送请求（请求100张照片的信息）后等待服务器响应会花费一些时间，于是就采用Python中的多线程，来同时发送多个请求。Threading是Python的一个标准模块，麻雀虽小，五脏俱全。单线程请求时，5000张照片的结果大约需要8分多钟，多线程请求时需要只需要40多秒即可。请求全部被相应后，即完成了所有照片信息的拼装工作：结果存储在类型为List的lphotos变量中。List是Python中使用最广泛的变量之一，灵活好用，一个顶5个~不过这里有个问题没有解决：在子线程的函数中，调用geoprocessor的addmessage()方法会报错（代码中注释掉的一句），非目前能力所能解决。
        题外话：上述多线程的工作也可以用并行计算来解决，这里有一篇出自Flyingis的关于Python中并行计算的文章，供大家学习。随着多核CPU的普及（ps：硅谷小公司推出100核CPU 性能是英特尔四倍），并行计算将会成为开发人员的新宠，微软.NET Framework 4.0中的Parallel Extention也说明了这一点。至于并行计算和多线程之间的区别，援引科班出身的同学胖胖鼠的话：

两者针对性不同。多线程技术是操作系统中的概念，是软件的概念。指进程中的一部分。以前并行计算主要更侧重硬件CPU，现在泛指同时使用多个计算资源解决问题。我感觉两者都用到了并发性并行性的概念。

个人感觉可简单理解为，并行计算是对细枝末节的多线程技术更优的抽象，与一个CPU里的多线程不同，它能利用多个CPU（不同服务器）资源，回归到网格计算的概念。下面继续主题。
        第二步：保存数据。之后我们便可以利用geoprocessor完成数据处理工作。首先通过调用Toolbox中的工具创建存储数据的FeatureClass，并调用工具添加相应的字段：

#create the 9.3 geoprocessor
gp=arcgisscripting.create(9.3)
gp.OverwriteOutput=True
myuserid=gp.getparameterastext(0)#"305320"
gp.Workspace=gp.getparameterastext(1)
myfile=os.path.basename(gp.getparameterastext(2))

#create 2 result featureclasses
try:
gp.addmessage("creating featureclasses...")
sr=gp.CreateObject("SpatialReference")
sr.Createfromfile(r"E:\Program Files\ArcGIS\Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj")
gp.Createfeatureclass(gp.workspace,myfile,"POINT","","","",sr)
gp.Addfield(myfile,"owner","TEXT","","",100)
gp.Addfield(myfile,"title","TEXT","","",100)
gp.Addfield(myfile,"photourl","TEXT","","",150)
#gp.Createfeatureclass(gp.workspace,xafile,"POINT",myfile,"","",sr)
del sr
except:
gp.adderror("creating featureclasses error")
raise Exception,"creating featureclasses error"

#push the photos into the featureclass
gp.addmessage("dealing my photos")
dealphotos(myuserid,myfile)

        最后两句便是整个功能的入口点了。dealphotos函数如下：

#function to push photos into featureclass
def dealphotos(set,fc):
photoscount=retrievephotos(set,0,1,"medium")
gp.addmessage("start to retrieving photos...")
getphotolist(set,0,photoscount,"medium")
gp.addmessage("retrieving photos completed")
cur=gp.InsertCursor(fc)
descfc=gp.Describe(fc)
print "begin insert features..."
gp.addmessage("begin insert features...")
for photo in lphotos:
row=cur.NewRow()
pnt=gp.CreateObject("POINT")
pnt.x=photo["longitude"]
pnt.y=photo["latitude"]
row.Setvalue(descfc.shapefieldname,pnt)
row.Setvalue("owner",photo["owner_name"])
row.Setvalue("title",photo["photo_title"])http://www.blogger.com/post-edit.g?blogID=21848155&postID=4114083803361112294#
row.Setvalue("photourl",photo["photo_file_url"])
cur.InsertRow(row)
print "insert features completed"
gp.addmessage("insert features completed")
del cur,descfc

        好啦，Python部分结束，添加到ArcMap中，便可成为一个Script Tool，与他人分享了。运行界面：

        运行结果：

        似乎不明显，那就利用空间分析中的Point Density工具抽象一下位置分布信息：

        也可以利用Model，将两部分连接起来，做成一个工具来用：

        有了地理数据，就可以继续在其他平台下进行处理和展示了：

最后总结一下标题：通过Python，不需要一个stand-alone的程序，就可以快速试验一个想法。
btw：附上最近google map tile的获取方法：

public override string GetTileUrl(int level, int row, int col)
{
int servernum = (col + 2 * row) % 4;
string sec2 = "&s="; // after &zoom=...
string secword = "Galileo";
int seclen = ((col * 3) + row) % 8; // corrected from int seclen = ( col + 3 * row ) % 8;
sec2 += secword.Substring(0, seclen);

//google maps map:
//http://mt1.google.com/vt/lyrs=m@114&hl=zh-CN&x=13&y=5&z=4&s=Gali
//http://mt1.google.com/vt/lyrs=m@114&hl=zh-CN&x=13&y=6&z=4&s=Galil
//http://mt1.google.com/vt/lyrs=m@114&hl=zh-CN&x=13&y=8&z=4&s=Galileo
string baseUrl = "http://mt" + servernum.ToString() + ".google.com/vt/lyrs=m@114&hl=zh-CN&x=";
string url = baseUrl + col.ToString() + "&y=" + row.ToString() + "&z=" + level.ToString() + sec2;
return url;

////google maps satallite
////http://khm1.google.com/kh/v=49&x=3&y=5&z=4&s=Galile
//string baseUrl = "http://khm" + servernum.ToString() + ".google.com/kh/v=49&x=";
//string url = baseUrl + col.ToString() + "&y=" + row.ToString() + "&z=" + level.ToString() + sec2;
//return url;
}


地理信息的另一种表达方式--Cartogram地图

        先来看一则旧闻：
英科学家绘制新世界人口地图 中印两国最突出
        下面是几张新闻里的地图：


        第一张图中可以看出，目前时间上人口最多的中国和印度被表达的很夸张，但这种图形上的夸张恰恰很好的表达了两国人口与其他国家人口数量的相对关系。这种地图被称为Cartogram：

The geometry or space of the map is distorted in order to convey the information of this alternate variable.

        为了要强烈表达地图中某种属性信息，而将图形进行一些扭曲，而扭曲的重要原则就是不改变原图形的拓扑关系。
        那么在ArcGIS Desktop中如何做出这种地图呢？请到ArcScripts网站下载这个免费的工具，安装好后它会出现在ToolBox中。它由密歇根大学的 Mark Newman和Michael Gastner开发完成。压缩包中已经包含了详细的使用说明和实例数据，下面是自带数据做出的一张2007年世界各国GDP总量地图：

        再次提醒，当你想突出地图上的某个属性信息而制作一张令人印象深刻的专题图时，记得试试这个工具。别忘了它适用矢量、栅格两种数据，而且可以做出经过两个或多个属性变量影响的cartogram地图。
        利用ArcGIS Desktop和这个工具，我们不是也作出了英国科学家才做出的地图吗：）

从arcgisscripting到arcpy——Geoprocessing在前进

        仔细查看9.3.1中本地的Desktop帮助会发现，没有任何先兆的情况下，An overview of writing geoprocessing scripts主题介绍了Arcpy模块，与在线帮助中该主题的内容完全不一样。而这个Arcpy就是9.4中arcgisscripting的升级版。文档编写人员在本地帮助中为9.4打了一个广告：）
        从win32com到arcgisscripting，geoprocessor实现了跨平台；从gp.create()到gp.create(9.3)，脚本中许多对象得到了质的升级。而arcpy正是arcgisscripting的继承者，它将为我们带来包括数据分析，数据转换，数据管理，及工作自动化在内的更加智能化的编程体验。在ESRI 2009 UC的视频中也能看到，ArcGIS Desktop 9.4中将Python环境更紧密的结合到了ArcMap中来。这也再次说明了Python作为scientific programming language for GIS的重要性。
        Geoprocessing是ArcGIS产品的三大核心（Geodatabase，Geovisualization，Geoprocessing）理念之一，也是AGS中gp服务的基石。Desktop之所以强大，在很大程度上取决于Geoprocessing Framework的支撑。若要深入掌握Geoprocessing，除了反复操练ArcToolbox，还得写的一手好的Python Scripts才行。