Distancia a una cresta en función del ángulo (Python)

Question

No tuve éxito buscando en Google durante un tiempo, pero es probable que no conozca la terminología correcta. Además, dime si ves otro intercambio más adecuado para esto.

Dado : Una matriz booleana de 2 dimensiones que almacena si un punto de la cuadrícula se considera una cresta de montaña (para la referencia estoy utilizando el algoritmo propuesto en Koka et al., Procedia Computer Science 4, (2011), 216-221 ).

Objetivo : Para un punto arbitrario de la cuadrícula, necesito un algoritmo que devuelva la distancia al punto más cercano que se considere una cresta (es decir, que sea un valor verdadero en la matriz de crestas) en función del ángulo (relativo al norte, por ejemplo). Debería permitir especificar un binning de los ángulos.

Hasta ahora sólo podía imaginar métodos de fuerza bruta, pero supongo que hay algo más inteligente.

Ejemplo : Tengo un modelo de elevación (ver las curvas de nivel en la parte izquierda) y el algoritmo de crestas determina que los píxeles negros son crestas. Necesito un retorno como el esbozado en la parte derecha, es decir, la distancia al punto más cercano de la cresta (delta_r) como función del ángulo (0-2pi = N-E-S-W). Por supuesto, habrá que hacer algunas excepciones para las direcciones en las que no hay cresta, como el NORTE en el ejemplo.

¿Tiene sugerencias sobre algoritmos eficaces?

EDIT, algunos enfoques de solución : Como referencia: Con el aporte de radouxju se me ocurrieron las siguientes soluciones tentativas.

FUNCIONES DE AYUDA :

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)

ENFOQUE 1 : Entonces se podría intentar definir un binning y encontrar puntos de cresta dentro de ese 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

Esto da lugar a problemas cuando la cresta no cercana es el ángulo sólido definido por el contenedor pero se encuentra una cresta "más lejana". Véanse los dos puntos de la parte inferior izquierda de la segunda imagen. Además, con este enfoque resultó no ser tan sencillo asignar "ningún valor" a las direcciones en las que no hay cresta.

ENFOQUE 2 : De ahí que se me ocurriera una especie de ventana móvil que identificara el punto de cresta más cercano al último (moviéndose en el sentido de las agujas del reloj). El código no se ha limpiado, algunas parcelas de depuración están sin comentar.

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_nuber)+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

Esto da el siguiente resultado:

Lo cual parece funcionar mejor, y ahora también se asignan correctamente los no valores.

Answer 1

La solución dependerá de la biblioteca que utilices (fiona, shapely, geopanda...). Mi algo sugerido está cerca de la fuerza bruta, pero no veo mucho más eficiente :

para cada punto, basándose en las coordenadas X e Y: - calcular la distancia de cada cresta de la montaña a su punto

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

• calcular el acimut de la línea que une el punto de la cresta y su punto

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

Una vez hecho esto, haz un bucle sobre el rumbo para seleccionar los puntos dentro de un rumbo determinado, y calcula la distancia mínima.

Answer 2

¿Has probado Script GDAL_proximity.py ? Esto se ejecuta desde la línea de comandos y toma los archivos como entradas.