Distance to a ridge as function of angle (Python)

Question

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.

Answer 1

Solution will depend on the library that you use (fiona, shapely, geopanda...). My suggested algo is close to brute force, but I don't see much more efficient :

for each point, based on X and y coordinates: - compute the distance of each mountain ridge to your point

def Distance(x1,y1,x2,y2):
        return ((x1-x2)^2+(y1-y2)^2)^0.5

• compute the azimuth of the line joining the ridge point and your point

def Azimuth(x1,y1,x2,y2):
        degBearing = math.degrees(math.atan2((x2 - x1),(y2 - y1)))
        if (degBearing < 0):
            degBearing += 360.0
        return degBearing

once this is done, loop on the bearing to select the points within a given bearing, and compute the minimum distance.

Answer 2

Have you tried GDAL_proximity.py script? This runs from the command line and takes files as inputs.