- 积分
- 56367
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-6-21
- 最后登录
- 1970-1-1

|

楼主 |
发表于 2018-9-13 23:18:06
|
显示全部楼层
可以用topology模块的 intersects 函数判断两个几何图形是否相交,参考此脚本:
 - import mipylib.geolib.topology as tp
- # Read typhoon data file
- fn = 'D:/Temp/ascii/CH2015BST.txt'
- tf = open(fn)
- lons = []
- lats = []
- prss = []
- for line in tf:
- #print line
- data = line.split()
- pn = int(data[2])
- for i in range(pn):
- line = tf.readline()
- data = line.split()
- lat = float(data[2])
- lats.append(lat * 0.1)
- lon = float(data[3])
- lons.append(lon * 0.1)
- t = data[0]
- prs = float(data[4])
- prss.append(prs)
- lons.append(nan)
- lats.append(nan)
- prss.append(nan)
- # Plot
- axesm()
- lworld = shaperead('D:/Temp/map/country1.shp')
- geoshow(lworld, facecolor=[200,200,200])
- layer = plotm(lons, lats, linewidth=1)
- xlim(100, 210)
- ylim(0, 60)
- title('Typhoon pathway')
- #Add a rectangle polygon
- lat = array([20, 20, 10, 10, 20])
- lon = array([180, 190, 190, 180, 180])
- g1 = geoshow(lat, lon, displaytype='polygon', color=[150,230,230,230], edgecolor='b', size=2)
- #Count the number of typhoons cross the rectangle
- n = 0
- for s in layer.shapes():
- if tp.intersects(s, g1):
- n += 1
- print 'Number: %i' % n
|
|