No success googling this for a while, but likely I just don't know correct terminology. Also, tell me if you see another exchange more suited for this.
Given: A boolean 2-dim array that stores if a point on the grid is considered a mountain ridge (for Reference I am using algorithm proposed in Koka et al., Procedia Computer Science 4, (2011), 216–221).
Goal: For an arbitrary point on the grid, I need an algorithm that returns the distance to the closest point that is considered a ridge (i.e. is a True value on the ridge-array) depending on the angle (relative to north e.g.). It should allow to specify a binning of the angles.
So far I could only imagine brute-force'ish step-methods, but I guess there is something cleverer.
Example: I have an elevation model (see contour lines on the left hand side) and the ridge alorithm deterimed the black pixels to be ridges. I need a return as sketched (!sry) on the right hand side, i.e. the distance to the closest ridge point (delta_r) as funciton of angle (0-2pi = N-E-S-W). Of course some exeption handling will be required for directions where there is no ridge i.e. NORTH in the example.
Do you have suggestions for efficient algorithms?
EDIT, some solution approaches: For reference: With radouxju's input I came up with the following tentative solutions.
HELPER FUNCTIONS:
def bearing(pole=[0,1], vec=[0,1]):
''' Calculate angle between pole and a vector.'''
ang1 = np.arctan2(*pole[::-1])
ang2 = np.arctan2(*vec[::-1])
return (ang1 - ang2) % (2 * np.pi)
def ridge_coords(east_grid,north_grid,ridges):
''' Find list of ridge coordinates.'''
ridge_indices = np.where(ridges==True)
ridges_east = east_grid[ridge_indices]
ridges_north = north_grid[ridge_indices]
ridge_coords = np.asarray((ridges_east,ridges_north)).T
return ridge_coords, ridge_indices
def dist2ridges(r_coords,location):
''' Calculate distance between an array of coords and a location.'''
return np.linalg.norm(r_coords-location,axis=1)
APPROACH 1: Then one could try to define a binning and find ridge points within that bin.
def radial_closest(bearings, dist, binning=np.linspace(0,2*np.pi,64)):
''' Find the radially closest ridge within an angle bin.'''
bin_idx = np.digitize(bearings,binning)
# print(np.hstack((bearings.reshape(-1,1),bin_idx.reshape(-1,1),dist.reshape(-1,1))))
all_data = np.ma.array(np.hstack((bearings.reshape((-1,1)),dist.reshape((-1,1)),bin_idx.reshape((-1,1)))),mask=False)
min_idx = np.ma.array(np.zeros_like(binning), mask=True)
min_dist = np.ones_like(binning)*np.nan
binning=binning+(binning[1]-binning[0])/2
for idx in range(len(binning)):
all_data.mask=False
all_data.mask=np.array([all_data[:,2]==idx]*3).T # true where bearings in bin
# all_data.mask=mask
if np.ma.is_masked(all_data):
all_data.mask=~all_data.mask
min_idx.mask[idx]=False
min_idx[idx]=np.argmin(all_data[:,1])
min_dist[idx]=np.amin(all_data[:,1])
min_idx=min_idx.astype(int)
return min_idx, min_dist, binning
def closest_ridges1(east_grid,north_grid,ridges_grid,location):
''' Put it all together, sort w.r.t. bearing and "format" return.'''
r_coords, r_indices = ridge_coords(east_grid,north_grid,ridges_grid)
dist = dist2ridges(r_coords,location)
bearings = np.array([bearing(vec=x-location) for x in r_coords])
min_idx, min_dist, binning = radial_closest(bearings, dist)
sorting = np.argsort(bearings[min_idx])
dist_interpolated = dist[min_idx][sorting]
close_ridge_indices = np.asarray(r_indices).T[min_idx][sorting]
close_ridge_indices = tuple(np.array(close_ridge_indices).T)
return binning, dist_interpolated, close_ridge_indices
This leads to problems when the no close ridge is the solid angle defined by the bin but a ridge "farther away" is found. See the two points on the lower left side of the second image. In addition with this approach it turned out to be not so simple to assign "no value" to directions where there is no ridge.
APPROACH 2: Hence I came up with sort of a moving window that identifieds the closest ridge point in vicinity of the last one (moving clockwise). Code not cleaned, some debuging plots are uncommented.
def closest_ridges2(east_grid,north_grid,ridges_grid,location):
r_coords, r_indices = ridge_coords(east_grid,north_grid,ridges_grid)
r_indices = np.asarray(r_indices).T
dist = dist2ridges(r_coords,location)
bearings = np.array([bearing(vec=x-location) for x in r_coords])
sorting=np.argsort(bearings)
close_ridges=[]
ridge_number=[]
all_data = np.hstack((bearings.reshape(-1,1),dist.reshape(-1,1),r_coords.reshape(-1,2),r_indices.reshape(-1,2)))[sorting,:]
all_data = np.hstack((all_data,np.arange(0,len(dist)).reshape(-1,1)))
all_data_copy = copy.deepcopy(all_data)
all_data = np.ma.array(all_data,mask=False)
min_dist_idx = all_data_copy[np.argmin(all_data[:,1]),6].astype(int)
close_ridges.append(min_dist_idx)
last_ridge = all_data_copy[close_ridges[-1],2:4]
next_in_all_data=close_ridges[-1] # index of
# print('closest ridge: {}'.format(close_ridges[-1]))
# print('coords: {}'.format(last_ridge))
#
radius=50
i=0
ridge_counter=0
while True:
# entering the loop the last_ridge element still is top
# print('loop number',i)
# print('ridge number', ridge_counter)
all_data.mask=False
# print('rolling up by: {}'.format(next_in_all_data))
all_data=np.roll(all_data,-next_in_all_data,axis=0)
all_data = np.delete(all_data, (0),axis=0)
dist_ = dist2ridges(all_data[:,2:4],last_ridge)
cond_dist = dist_> radius
bearings_ = np.array([bearing(pole=last_ridge-location,vec=vec-location) for vec in all_data[:,2:4]])
cond_bearings = bearings_ >= np.pi
cond = np.logical_or(cond_bearings,cond_dist)
ridge_number.append(ridge_counter)
if np.all(cond_dist):
ridge_counter+=1
#
all_data.mask =np.array([cond]*7).T
next_in_all_data = np.argmin(all_data[:,1])
all_data.mask=False
min_dist_idx = all_data[next_in_all_data,6] # copy of external index
close_ridges.append(min_dist_idx)
last_ridge = all_data[next_in_all_data,2:4]
# print('next point in all data: {}'.format(next_in_all_data))
# if next_in_all_data == 0: # keep
# pass # keep
# elif next_in_all_data > 0: # keep
all_data=np.delete(all_data, slice(0,next_in_all_data),axis=0) #apparently the deletion sets the mask to false (and deletes masked as other rows)
next_in_all_data=0
# fig,ax=plt.subplots()
# fig.show
#
# ax.scatter(all_data[:,2],all_data[:,3],marker='^',c='g')
# ax.scatter(all_data_copy[np.array(close_ridges,dtype=int),2],all_data_copy[np.array(close_ridges,dtype=int),3],marker='x',c='r')
# ax.scatter(location[0],location[1],marker='X',c='k')
# ax.add_patch(plt.Circle((location[0], location[1]), radius, color='r', alpha=0.2,linestyle='--'))
# ax.add_patch(plt.Circle((last_ridge[0], last_ridge[1]), radius, color='r', alpha=0.2,linestyle='--'))
# ax.text(x=0.5,y=0.05,s=str(i)+' '+str(np.all(cond_dist))+' '+str(ridge_counter),transform=ax.transAxes)
if all_data.shape[0]==1:
# print('break')
break
i+=1
ridges_list=[[] for x in range(max(ridge_number)+1)]
bearing_list=[[] for x in range(max(ridge_number)+1)]
distance_list=[[] for x in range(max(ridge_number)+1)]
for x,y in zip(close_ridges,ridge_number):
ridges_list[y].append(int(x))
bearing_list[y].append(all_data_copy[int(x),0])
distance_list[y].append(all_data_copy[int(x),1])
print(bearing_list)
binning = np.linspace(0,2*np.pi,256)
dist_interpolated=np.zeros_like(binning)
for b,d in zip(bearing_list,distance_list):
if b[0]<b[-1]:
inter_ = np.interp(binning,b,d,left=0,right=0)
if b[0]>b[-1]:
inter_ = np.interp(binning,b,d,period=2*np.pi)
inter_[np.logical_and(binning>b[-1],binning<b[0])]=0
dist_interpolated=dist_interpolated+inter_
dist_interpolated[dist_interpolated == 0] = np.nan
# print(np.hstack((binning.reshape(-1,1),dist_interpolated.reshape(-1,1))))
# fig,ax=plt.subplots()
# for i,plotlist in enumerate(ridges_list):
# ax.scatter(all_data_copy[plotlist,0],all_data_copy[plotlist,1],marker='x',label='ridge nr {}'.format(i))
#
# ax.plot(binning,dist_interpolated,c='k',linestyle='--',linewidth=4)
# for b,d in zip(bearing_list,distance_list):
# ax.plot(b,np.array(d),marker='o')
close_ridge_indices = tuple(np.array(all_data_copy[np.array(close_ridges,dtype=int),4:6],dtype=int).T)
return binning, dist_interpolated, close_ridge_indices
This yields the following result:
Which appears to work better, and now also non-values are properly assigned.