|
USGS官方網(wǎng)站每天都會(huì)實(shí)時(shí)更新全世界的地震信息，包含地震發(fā)生的地點(diǎn)，坐標(biāo)，震級(jí)，震中距離地表的距離等等，坐標(biāo)系采用的是地心坐標(biāo)系WGS84，如何將這些實(shí)時(shí)的信息采集到自己的系統(tǒng)之中，用于進(jìn)一步的科學(xué)計(jì)算和空間分析，需要借助一些解析數(shù)據(jù)和空間計(jì)算的方法。 GIS主流的應(yīng)用策略之一，是融合共享，也是技術(shù)發(fā)展的整體需求，接下來(lái)的內(nèi)容并不關(guān)心如何共享，而是從共享的最基礎(chǔ)層面-數(shù)據(jù)層面，來(lái)解析USGS網(wǎng)站頁(yè)面中的地震空間數(shù)據(jù)。 地震信息源的網(wǎng)址：http://quake.wr./recenteqs/Quakes/quakes0.htm 抽取該網(wǎng)站上的地震數(shù)據(jù)，可以使用主流的.Net，Java，開(kāi)發(fā)一個(gè)小程序，相對(duì)獨(dú)立，但是更多的空間分析和應(yīng)用是直接基于GIS桌面平臺(tái)進(jìn)行的，如同在 MatLab平臺(tái)上實(shí)現(xiàn)了一套數(shù)學(xué)分析的思路，需要引入更多的外部資源充實(shí)其中的變量和數(shù)組，因此，我們可以直接在GIS桌面平臺(tái)ArcMap中直接用 Python去抽取數(shù)據(jù)，不同于普通的變量和數(shù)組，空間數(shù)據(jù)的引入還需要考慮坐標(biāo)轉(zhuǎn)換，符號(hào)設(shè)置等相關(guān)信息，基于AO的Python都可以幫助咱們?nèi)ヒ灰粚?shí)現(xiàn)。 urllib模塊是標(biāo)準(zhǔn)Python庫(kù)的一部分，方便提取最原始的地震數(shù)據(jù)：  for l in q.readlines(): if l.find("<STRONG>") == 0: l = l[8:] l = l.replace("FONT COLOR", "FC") if l.find('<A HREF="/recenteqs/Maps') == 0: quakeI = l.split() magnitude = float(quakeI[2]) x = -float(quakeI[7][:-1]) y = float(quakeI[6][:-1]) point = arcpy.Point(x,y) feature = cur.newRow() feature.shape = point feature.setValue("magnitude", magnitude) cur.insertRow(feature) 這樣所有的屬性信息，包括x/y坐標(biāo)數(shù)據(jù)都已經(jīng)獲取。熟悉AO的朋友肯定非常了解"cur.insertRow(feature)"的含義和開(kāi)發(fā)過(guò)程了。不熟悉的請(qǐng)看后面解釋吧：
SR = arcpy.SpatialReference(r"C:\workspace\demo3\WGS 1984.prj") #創(chuàng)建一個(gè)FeatureClass arcpy.CreateFeatureclass_management(os.path.dirname(tmpFC), os.path.basename(tmpFC), "POINT", spatial_reference = SR) #增加一個(gè)屬性字段，代表地震震級(jí) arcpy.AddField_management(tmpFC, "Magnitude", "DOUBLE") #獲取該FeatureClass的插入游標(biāo) cur = arcpy.InsertCursor(tmpFC) 所以"cur.insertRow(feature)"就是將所有網(wǎng)站上獲取的每一個(gè)元組信息都添加到tmpFC臨時(shí)FeatureClass之中了。
現(xiàn)在數(shù)據(jù)已經(jīng)獲取了，按照傳統(tǒng)的解析數(shù)據(jù)方法，咱們?nèi)蝿?wù)也就完成了，但是對(duì)于GIS應(yīng)用來(lái)說(shuō)，需要將這些數(shù)據(jù)顯示到基礎(chǔ)地圖上，這里面就需要思考兩個(gè)問(wèn)題： 1.已獲取數(shù)據(jù)的坐標(biāo)系和基礎(chǔ)地圖數(shù)據(jù)是否相同？不相同則需要坐標(biāo)轉(zhuǎn)換。 2.如何符號(hào)化顯示？ 假如我們基礎(chǔ)地圖的坐標(biāo)是North_America_Albers_Equal_Area_Conic.prj，坐標(biāo)轉(zhuǎn)換可以通過(guò)以下兩行代碼完成： SR = arcpy.SpatialReference(r"C:\workspace\demo3\North_America_Albers_Equal_Area_Conic.prj") arcpy.Project_management(tmpFC, outFC, SR, "NAD_1983_To_WGS_1984_1") 符號(hào)化可以自定義，也可以參考已有圖層的樣式，如： arcpy.ApplySymbologyFromLayer_management(os.path.splitext(os.path.basename(outFC))[0], r"c:\workspace\demo3\earthquake.lyr") 這樣咱們就完成了解析的工作，可以在ArcMap中基于這些數(shù)據(jù)進(jìn)行進(jìn)一步的分析。ArcMap 9.3需要在單獨(dú)的IDLE環(huán)境中寫(xiě)開(kāi)發(fā)腳本，將結(jié)果手工添加到ArcMap平臺(tái)軟件中，ArcMap 9.4則直接整合了Python運(yùn)行環(huán)境，開(kāi)發(fā)過(guò)程中加入了動(dòng)態(tài)提示，實(shí)時(shí)幫助，應(yīng)用交互等等，在科學(xué)計(jì)算和空間分析中非常方便。
插圖： 原始軟件界面 導(dǎo)入或直接編寫(xiě)Python代碼 右側(cè)幫助信息 |
|
|
來(lái)自： LibraryPKU > 《WebGIS》
