我有两个纬度和经度列表,分别代表构成地质断层的一堆点。我正在尝试制作一个通过底图(通过底图)将这些故障显示为线的地图,但是,使用pyplot.plot函数,这些点将连接相距数英里的故障,我不希望这样。
到目前为止,我一直在做的只是使这些点变大,以致看起来好像是一堆线,但是当我放大时,显然是点而不是线。
知道组成单个断层的这些点之间相隔50米可能也很有用。最初,我有一个代码,使用大圆工具确定了这些点之间的距离,并创建了一个新的距离列表:
f_distance = []
c = 0
for i,j in zip(fault_lat,fault_lon):
while c < (len(fault_lon)-1):
location1 = (fault_lat[c],fault_lon[c])
location2 = (fault_lat[c+1],fault_lon[c+1])
distance = great_circle(location1,location2).meters
f_distance.append(distance)
c+=1
然后,如果该距离小于或等于50,它将绘制它们并(希望)绘制它们之间的线,否则将跳过它。
c = 0
for i,j,k in zip(map_q_fault_lon,map_q_fault_lat,q_fault_distance):
if k <= 50:
map.plot(i,j,linestyle = '-',linewidth = 1.00,color = 'black')
c+=1
但是可惜的是它给了我一个空白的屏幕。无论如何,我是否可以对此进行重新加工,使其能够满足我的要求?
对于约50m的距离,您实际上不需要很大的圆弧距离就能获得良好的近似值。我假设fault_lat
和fault_lon
是弧度的numpy数组。如果没有,请转换它们:
fault_lat = np.radians(fault_lat)
fault_lon = np.radians(fault_lon)
现在,您可以使用简单的公式来计算距离,该公式涉及地球的半径和适当的纬度缩放比例:
r = 6378000
dist = r * np.sqrt(np.diff(fault_lat)**2 + (np.diff(fault_lon) * np.cos(0.5 * (fault_lat[1:] + fault_lat[:-1])))**2)
这将创建连续点之间(近似)距离的索引。您可以在距离超过某个阈值(例如50m)的位置拆分结果:
threshold = 50
indices = np.flatnonzero(dist > threshold) + 1
lat_segments = np.split(fault_lat, indices)
lon_segments = np.split(fault_lon, indices)
您可以使用以下方法绘制结果:
for lat, lon in zip(lat_segments, lon_segments):
map.plot(lon, lat)