USGS官方网站每天都会实时更新全世界的地震信息,包含地震发生的地点,坐标,震级,震中距离地表的距离等等,坐标系采用的是地心坐标系WGS84,如何将这些实时的信息采集到自己的系统之中,用于进一步的科学计算和空间分析,需要借助一些解析数据和空间计算的方法。
GIS主流的应用策略之一,是融合共享,也是技术发展的整体需求,接下来的内容并不关心如何共享,而是从共享的最基础层面-数据层面,来解析USGS网站页面中的地震空间数据。
地震信息源的网址:http://quake.wr.usgs.gov/recenteqs/Quakes/quakes0.htm
抽取该网站上的地震数据,可以使用主流的.Net,Java,开发一个小程序,相对独立,但是更多的空间分析和应用是直接基于GIS桌面平台进行的,如同在 MatLab平台上实现了一套数学分析的思路,需要引入更多的外部资源充实其中的变量和数组,因此,我们可以直接在GIS桌面平台ArcMap中直接用 Python去抽取数据,不同于普通的变量和数组,空间数据的引入还需要考虑坐标转换,符号设置等相关信息,基于AO的Python都可以帮助咱们去一一实现。
urllib模块是标准Python库的一部分,方便提取最原始的地震数据:
q
=
urllib.urlopen(r
'
http://quake.wr.usgs.gov/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)
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坐标数据都已经获取。
熟悉AO的朋友肯定非常了解"cur.insertRow(feature)"的含义和开发过程了。不熟悉的请看后面解释吧:
#
定义WGS84坐标参考系
SR = arcpy.SpatialReference(r " C:\workspace\demo3\WGS 1984.prj " )
# 创建一个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)
SR = arcpy.SpatialReference(r " C:\workspace\demo3\WGS 1984.prj " )
# 创建一个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)"就是将所有网站上获取的每一个元组信息都添加到tmpFC临时FeatureClass之中了。
现在数据已经获取了,按照传统的解析数据方法,咱们任务也就完成了,但是对于GIS应用来说,需要将这些数据显示到基础地图上,这里面就需要思考两个问题:
1.已获取数据的坐标系和基础地图数据是否相同?不相同则需要坐标转换。
2.如何符号化显示?
假如我们基础地图的坐标是North_America_Albers_Equal_Area_Conic.prj,坐标转换可以通过以下两行代码完成:
现在数据已经获取了,按照传统的解析数据方法,咱们任务也就完成了,但是对于GIS应用来说,需要将这些数据显示到基础地图上,这里面就需要思考两个问题:
1.已获取数据的坐标系和基础地图数据是否相同?不相同则需要坐标转换。
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.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 " )
r " c:\workspace\demo3\earthquake.lyr " )
这样咱们就完成了解析的工作,可以在ArcMap中基于这些数据进行进一步的分析。ArcMap 9.3需要在单独的IDLE环境中写开发脚本,将结果手工添加到ArcMap平台软件中,ArcMap 9.4则直接整合了Python运行环境,开发过程中加入了动态提示,实时帮助,应用交互等等,在科学计算和空间分析中非常方便。
插图:
原始软件界面
自动提示
插图:
原始软件界面
最终结果