解析USGS網(wǎng)站頁面中的地震空間數(shù)據(jù)

LibraryPKU 2019-10-28

展開全文

作者：Flyingis
本文歡迎友情轉載，但請注明作者及原文鏈接，嚴禁用于商業(yè)目的！

    USGS官方網(wǎng)站每天都會實時更新全世界的地震信息，包含地震發(fā)生的地點，坐標，震級，震中距離地表的距離等等，坐標系采用的是地心坐標系WGS84，如何將這些實時的信息采集到自己的系統(tǒng)之中，用于進一步的科學計算和空間分析，需要借助一些解析數(shù)據(jù)和空間計算的方法。

    GIS主流的應用策略之一，是融合共享，也是技術發(fā)展的整體需求，接下來的內容并不關心如何共享，而是從共享的最基礎層面-數(shù)據(jù)層面，來解析USGS網(wǎng)站頁面中的地震空間數(shù)據(jù)。

    地震信息源的網(wǎng)址：http://quake.wr./recenteqs/Quakes/quakes0.htm

    抽取該網(wǎng)站上的地震數(shù)據(jù)，可以使用主流的.Net，Java，開發(fā)一個小程序，相對獨立，但是更多的空間分析和應用是直接基于GIS桌面平臺進行的，如同在 MatLab平臺上實現(xiàn)了一套數(shù)學分析的思路，需要引入更多的外部資源充實其中的變量和數(shù)組，因此，我們可以直接在GIS桌面平臺ArcMap中直接用 Python去抽取數(shù)據(jù)，不同于普通的變量和數(shù)組，空間數(shù)據(jù)的引入還需要考慮坐標轉換，符號設置等相關信息，基于AO的Python都可以幫助咱們去一一實現(xiàn)。

    urllib模塊是標準Python庫的一部分，方便提取最原始的地震數(shù)據(jù)：

q = urllib.urlopen(r'http://quake.wr./recenteqs/Quakes/quakes0.htm')
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坐標數(shù)據(jù)都已經(jīng)獲取。熟悉AO的朋友肯定非常了解"cur.insertRow(feature)"的含義和開發(fā)過程了。不熟悉的請看后面解釋吧：

    #定義WGS84坐標參考系
    SR = arcpy.SpatialReference(r"C:\workspace\demo3\WGS 1984.prj")
    #創(chuàng)建一個FeatureClass
    arcpy.CreateFeatureclass_management(os.path.dirname(tmpFC),
                                        os.path.basename(tmpFC),
                                        "POINT",
                                        spatial_reference = SR)
    #增加一個屬性字段，代表地震震級
    arcpy.AddField_management(tmpFC, "Magnitude", "DOUBLE")
    #獲取該FeatureClass的插入游標
    cur = arcpy.InsertCursor(tmpFC)

    所以"cur.insertRow(feature)"就是將所有網(wǎng)站上獲取的每一個元組信息都添加到tmpFC臨時FeatureClass之中了。

    現(xiàn)在數(shù)據(jù)已經(jīng)獲取了，按照傳統(tǒng)的解析數(shù)據(jù)方法，咱們任務也就完成了，但是對于GIS應用來說，需要將這些數(shù)據(jù)顯示到基礎地圖上，這里面就需要思考兩個問題：

    1.已獲取數(shù)據(jù)的坐標系和基礎地圖數(shù)據(jù)是否相同？不相同則需要坐標轉換。
    2.如何符號化顯示？

    假如我們基礎地圖的坐標是North_America_Albers_Equal_Area_Conic.prj，坐標轉換可以通過以下兩行代碼完成：

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

符號化可以自定義，也可以參考已有圖層的樣式，如：

arcpy.ApplySymbologyFromLayer_management(os.path.splitext(os.path.basename(outFC))[0],
r"c:\workspace\demo3\earthquake.lyr")

    這樣咱們就完成了解析的工作，可以在ArcMap中基于這些數(shù)據(jù)進行進一步的分析。ArcMap 9.3需要在單獨的IDLE環(huán)境中寫開發(fā)腳本，將結果手工添加到ArcMap平臺軟件中，ArcMap 9.4則直接整合了Python運行環(huán)境，開發(fā)過程中加入了動態(tài)提示，實時幫助，應用交互等等，在科學計算和空間分析中非常方便。

    插圖：

    原始軟件界面